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Abstract 
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s modeling was developed under this grant monitored by Dr. Peter 
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period is November 14, 1986 to January 31, 1993. 
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1 Introduction 

In the study of high speed flows at high altitudes, such as that encoun- 
tered by re-entry spacecrafts, the interaction of chemical reactions and other 
non-equilibrium processes in the flow field with the gas dynamics is cru- 
cial. Generally speaking, problems of this level of complexity must resort 
to numerical methods for solutions, using sophisticated computational fluid 
dynamics (CFD) codes. The difficulties introduced by reacting gas dynamics 
can be classified into three distinct headings: 

I. The usually inadequate knowledge of the reaction rate coefficients in the 

non-equilibrium reaction system, and * 

II. The vastly larger number of unknowns involved in the computation and 

the expected stiffnesd^ of the equations, and 

III. The interpretation of the detailed reacting CFD numerical results. 

It is not uncommon to have to deal with reaction systems involving tens of 
reactants which participate in hundreds of elementary reactions (using spec- 
ulative rate coefficients) which have vastly disparate time scales. Assuming 
that sufficient computing resources are available, that the most up-to-date re- 
action rate coefficients have been used, we must still address item #111 above: 
how does one extract physical understanding from the massive amount of nu- 
merical results generated by a validated reacting CFD code? 

Traditionally, an important component of a theory is its formulation. By 
limiting its validity on some specific domain of the parameter space, a tradi- 
tional theoretician “neglects” the unessential complications and retains only 
the barest mechanisms to arrive at the simplified model which nevertheless 
fully mimics the behaviors of the complete full-model. In most cases, it is 
the simplified model itself, rather than the detailed (analytical or numerical) 
solution of the simplified model, that offers the most physical insights to the 
problem under investigation. This traditional approach is distinctly different 
from the modern CFD approach: retain all complications which can be over- 
come by massive computing resources; just generate the numerical solutions, 
then worry about what physical insights can be extracted from the solutions 
later. 

The research performed under this grant accepts the premise that react- 
ing flows of practical interest in the future will in general be too complex or 



Final Report for NASA Grant NAG-1-726 


3 


“untractable” for traditional analytical developments. The power of modern 
computers must be exploited. However, instead of focussing solely on the 
construction of numerical solutions of full-model equations, attention is also 
directed to the “derivation” of the simplified model from the given full- model. 
In other words, the present research aims to utilize computations to do tasks 
which have traditionally been done by skilled theoreticians: to reduce an 
originally complex full-model system into an approximate but otherwise e- 
quivalent simplified model system. The tacit assumption is that once the 
appropriate simplified model is derived, the interpretation of the detailed 
numerical reacting CFD numerical results will become much easier. 

The approach of the research is called computational singular perturbation 
(CSP) 1 [2] which is developed under the present grant. Dr. D. A. Goussis, 
currently assistant professor at University of Patra in Greece, was a research 
staff member and made significant contributions to the work reported here. 


2 Goals of the Research 

The traditional approach to derive simplified kinetics models uses the so- 
called steady-state and partial- equilibrium approximations™ U 5] ’ 1 u ],t l The 
steady-state approximation is applied to certain radicals or reaction inter- 
mediaries,” and the partial-equilibrium approximation is applied to certain 
“fast” reactions. For relatively simple and therefore “tractable” reaction sys- 
tems, these approximations are applied guided by experience and intuition. 
The results of such analytical studies is a “reduced” reaction system. The 
concentration of the radicals can be found from a set of algebraic equations 
of state derived from the approximations, enabling the reduction of the num- 
ber of unknowns. The simplified kinetics equations derived is no longer stiff. 
Most importantly, much insight can be gain from inspecting the equations 
of state and the simplified kinetics equations. Usually, once an analytical 
simplified model is obtained, the following list of questions can be answered 
intelligently by inspection of the results. 

1. How many algebraic "equations of states” are available and what are 
they? 

2. How many "conservation laws” are available and what are they? 



Final Report for NASA Grant NAG-1-726 


4 


3. What reactions are so fast that their reaction rate coefficients need not 
be known accurately so long as they are sufficiently fast? 

4. What reactions are so slow that they can be ignored altogether? 

5. What reactions are primarily responsible for the current behavior of 
the solution? 

6. What is the effective stoichiometry of the reaction system? 

7. Which reactants can be solved for in terms of the others from the 
available equations of state? 

8. What reactants can be summarily excluded from the reaction system? 

9. What is the response of the reaction system to some perturbation? 

Given a massively complex chemical reaction system involving N un- 
knowns and R elementary processes where N and R are large numbers, the 
traditional approach is not viable. The goal of the present research is to 
generate the simplified model and the answers to the above questions using 
numerical computations. 


3 Example 


Consider a simple reaction system with three species A, B and C and two 
elementary reactions: 


Reaction #1 


A + A^B 


(3.1) 


Reaction #2 


A- B + C 


(3.2) 


The two unknowns can be considered the components of a column vector 
y = [A,B] t . The stoichiometric vectors s r and reaction rates F r of the 
elementary reactions are expressed as follows: 

s, = I— 2,1, Of, F' = kM' -K*B), (3.3) 

Sj = [-1,1, If, F 2 = h(A - K 2 BC), (3.4) 
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where k u k 2 , K x and K 2 are rate coefficients obtained from up-to-date and 
reliable databases. The governing equations can be written as: 

% = s,F' + SjF J . (3.5) 

at 

Now suppose it is known that reaction #2 is very fast and becomes exhausted 
very quickly. How does one take advantage of this observation to derive the 
simplified model equations? 

Obviously, when reaction #2 is exhausted, one would have: 

F 2 = k 2 (A - K 2 BC) « 0 (3.6) 

which is an algebraic equation of state relating A, B and C. However, it is 
easy to verify that if one substitutes (3.6) back into (3.5), one would obtain 
the wrong answer for the simplified model equations. The difficulty of this 
simple-minded but incorrect derivation is quite subtle but is widely known. 
Consequently, special ad hoc procedures are recommended by various au- 
thors on how to proceed once the fast reactions of a reaction system have 
been identified. A common feature of these ad hoc procedures is that the 
unknowns are divided into two groups: those which can be eliminated from 
the differential equation system by the use of the derived equations of state, 
and those which cannot. The first group is usually called radicals or inter- 
mediaries. For the example here, a knowledgeable kineticist will correctly 
conclude that either A or B can be treated as radicals, while C cannot. 

Normally, a chemical reaction system respects the Law of Conservation of 
atomic elements. While such conservation laws are obviously automatically 
satisfied by a correctly formulated model, their existence is usually identified 
on physical grounds rather than on mathematical grounds. For the present 
example which is hypothetical, the physical identification of such laws become 
impossible. However, it can readily be verified mathematically that A + 2B- 
C is the only conserved scalar for this reaction system. 

We can represent (3.5) by the following alternative representation: 

^ « sj.F 1 ' + s r F*' (3-7) 

at 

where the primed reactions are some linear combination of the original ele- 
mentary reactions: 
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Reaction #1’ 


Si' = Sj - cs 2 , 


F 1 ' = F\ 


(3.8) 


Reaction #2’ 


s 2 > — s 2 , 


F 2 ' = F 2 + cF\ 


(3.9) 


where c is an arbitray constant. Note that this alternative representation is 
exact, and that the primed reactions do not, in general, have any “physical 
meanings.” 

The conventional steady-state approximation in essence recommends that 
c be chosen such that the radical component of the effective stoichiometric 
vector of the primed reaction driven solely by the slower original reactions 
be zero. Hence for our example we would choose c = 2 if A is identified 
as the radical (i.e. s r = [0,-1, -2] r ), while c = 1 if B is so identified 
(i.e. s 2 < = [-1,0, -l] r ). In addition, the resulting differential equation for 
the radical is discarded and is replaced by the equation of state obtained by 
setting the reaction rate of the fast primed reaction to zero. For our example, 
this equation of state is: 

F 2 ' = F 2 + cF x = 0. (3.10) 


Using this derived equation of state and the conservation of A + 2B - C, 
a simplified model of the original reaction system is derived, consisting of a 
single differential equation. As can be seen from the above derivation, the re- 
sults obtained depends on the identification of the radical. In addition to the 
above described steady-state approximation, a different set of slightly more 
complicated ad hoc procedures are recommended for the partial-equilibrium 
approximation. It can be shown that each of these procedures generates a 
different analytical result, and the validity and accuracy of these results can 
only be assessed empirically. 

In general, it requires experience, intuition and skill in order to make the 
correct derivation of the simplified models using the conventional approach. 


4 The Basic Ideas of CSP 

Consider a reaction system consisting of N unknowns (species, temperature 
and density) and R elementary reactions. Limiting ourselves to spatially 
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homogeneous problems, the governing system of ODE s can be expressed as 
follows: , 

Tt = g(y) (4 ' I) 

where yand gare iV-dimensional column vectors representing the unknowns 
and their time rate of change, respectively. In general, the full-model g can 
be expressed as the sum of R terms each representing a physically meaningful 
elementary process: 

g = I> r F r (y) ( 4 - 2 ) 

r= 1 

where s r and F r (y) are the stoichiometric vector and the reaction rate of the 
r-th elementary process, respectively. In general, N and R are different, and 
for most full-model reaction systems R is usually larger than ?V. Equation 
(4.2) is called the physical representation of g, because each additive term 
can be associated with some physical meaning by the investigator. 

The basic idea of CSP is to expresses g in terms of a set of N (yet to be de- 
termined) linearly independent column basis vectors aj(t), ( i = 1,2, .. . , N). 

g = X>/‘ (43) 

t = l 

where /’ is the “amplitude” of the z'-th mode given by: 

r = b*.g, *■ = 1,2,..., W, (4.4) 

and b'(t), (i=l, 2, . . . , N) are row vectors which are orthonormal to aj(t): 

b i -a j = <5}., ij = 1,2,..., W. (4.5) 

Obviously, the column basis vectors aj’s are some linear combinations of the 

physically meaningful stoichiometric vectors s r ’s, and the mode amplitudes /' 
are some consistent linearly combinations of the original elementary reaction 
rates F r, s. The question asked by CSP is: how should the set of basis vectors 
aj be chosen such the “fast” modes can be neglected from (4.3) after they 
are exhausted? 

To answer this question in general, we need to know how the mode ampli- 
tudes /* evolve with time. Differentiating /* with respect to time, we obtain, 
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with the assistance of (4.5), the following equation: 

( 4 - 6 > 

j = 1 

where 

A‘ = (—• + b‘ • J) • aj. i,j = l,2,...,N. (4.7) 

If the N x N real matrix A’ were diagonal, the N modes would be completely 
decoupled; the amplitude of each mode would evolved according to its own 
time scale which can be identified with the reciprocal of the diagonal element. 
These modes can be ordered in ascending order of the magnitude of the time 
scales. Assuming that all modes are the decaying type, the amplitudes of the 
fast modes would decay toward zeros. At any moment of time, the simplified 
chemical kinetics model can be obtained by simply neglecting the fast modes 
which have become exhausted. 

But in general, the matrix A .*■ is not diagonal. The mathematical theory of 
CSP is focussed on finding basis vectors which minimizes the the off-diagonal 
terms of this matrix. These CSP-derived basis vectors are used to compute 
additional CSP data which can explicitly provide physical insights about the 
reaction system under study. 


5 Accomplishments 

The CSP theory is developed incrementally throughout the grant period. 
Appendix I is a reprint of the most recently published paper 101 and contains 
the most complete and up-to-date summary of the theoretical developments. 

We have mainly focussed on homogeneous reaction systems so that ef- 
fects of diffusion are not included. A paper 1101 was presented at the 1992 
24th International Combustion Symposium at Sydney on the oxidation of 
methanol. This full model reaction system consisted of 30 chemical species 
plus temperature as unknowns and involved 173 reversible reactions. The 
CSP-derived minimum reaction system consisted only 15 species plus tem- 
perature. A copy of this published manuscript is included here as Appendix 

II. 

A paper 1111 was presented at the 1992 APS Fluid Dynamics Division 
Meeting dealing specifically with the effects of transverse diffusion. A full 
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length paper of this work is in preparation, and a copy of the unfinished 
manuscript is included here as Appendix III. 

The CSP theory involves moderately abstract concepts in linear algebra 
and higher mathematics. A common comment on CSP from chemical kineti- 
cists is that the mathematics involved is too difficult. A full length paper 
which tries to explain CSP using simple examples and actual numbers is 
in preparation, and a copy of the unfinished manuscript is included here as 
Appendix IV. 

The following is a list of the major scientific accomplishments performed 
under this grant: 

Computer Code CSP91 A computer code, named CSP91, was written 
by Dr. D. A. Goussis and contained a CSP Fortran subroutine library 
which can be called by reacting flow codes. This code is compatible 
with Chemkin 11 '* 1 , and uses Chemkin database format and data files. 
Only minimum documentation is provided via a number of readme files. 

A copy of this code has been delivered to Dr. Peter Gnoffo of NASA 
Langley. 

Clarification of the Conventional Approaches The reliance on experi- 
ence and intuition in the conventional approaches to simplified kinetics 
modelling is completely removed. In essence, CSP considers the con- 
ventional approaches as guessing at the correct fast basis vectors. By 
showing that the ideal basis vectors should diagonalize A*-, it sets the 
goal for the ensuing theoretical work. It shows explicitly that the local 
eigen- vectors of the Jacobian of g will not diagonalize Aj J ^ 

The Refinement Procedure Using all N Basis Vectors The new idea 
of CSP is not to look for basis vectors which will precisely diagonalize 
A* . It proposes an iteration process via a refinement procedure. One 
may begin with a reasonable guess, and the refinement procedure will 
generate a better iterant. The eigen-vectors of the Jacobian of g is 
always available to provide an excellent first guess. 

The Refinement Procedure Using Only Fast Basis Vectors In all the 
published papers, the refinement procedure requires that all N basis 
vectors participate in the algorithm even though theoretically only the 
fast basis vectors are of interest. In the unpublished paper presented 
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in Appendix IV, a simpler and much cleaner refinement procedure is 
introduced in which only the fast basis vectors are involved. 

The CSP-derived Equations of State When the amplitudes of the de- 
caying fast modes are sufficiently small in some pragmatic sense, they 
become equations of states. One no longer needs to apply the steady- 
state approximation or the partial-equilibrium approximation to derive 
them — such equations of state can now be computationally derived rou- 
tinely. 

The CSP-derived Conservation Laws The concept of conservation laws 
are generalized. In any given time epoch, reaction modes which are 
too slow to be included in an approximate theory can be exploited 
to yield locally valid conserved scalars. In other words, conservation 
laws can now also be computationally derived, including the familiar 
conservation of atomic elements. 

The Radical Pointer Once CSP-derived basis vectors are available, a rad- 
ical pointer can be constructed which identifies which chemical species 
can be algebraically solved for in terms of the others using the available 
equations of state. One no longer needs to guess or argue which chem- 
ical species should be eliminated from the differential equation system. 
The radical pointer specifically identify the “wrong” species for this 
purpose. 

The Fast Reaction Pointer In general, equations of state must be used 
with great care. While they are doubtlessly valid approximate equa- 
tions of state, they cannot be used blindly to manipulate the differential 
equations. For example, one may not used the M available equations 
of state to eliminate any M original elementary reaction rates from 
the original system of differential equations. Only the “fast” original 
elementary reaction rates can be eliminated, and wrong answers would 
be obtained if the wrong reactions were eliminated. The theory of CSP 
provides a fast reaction pointer , allowing the fast reactions be identified 
computationally. 

Minimum Reaction System The code CSP91 written by Dr. D. A. Gous- 
sis can computationally generate from the original full-model reaction 
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system a minimum reaction system which contains a much smaller num- 
ber of species (and includes any species of interest specifically requested 
by the investigator). 

Effects on Diffusion Terms When diffusion is included in a reacting flow 
problem, the coupling between fast reactions and the diffusion process 
of the participants is strong^ 141 . In general, one must find the new 
form of the diffusion term, and in addition, one must also derived the 
new boudnary conditions consistent with the exhaustion of fast modes. 
Significant progress has been made on these issues, and the current 
status of this research is reported in Appendix III. 

CSP and The RNG Theory of Turbulence While our interest in de- 
veloping CSP was prompted by reacting gas flow problems, the central 
mathematical issue of the theory is how to deal with non-linear prob- 
lems with vastly disparate time scales. The chemical kinetics problems 
are characterized by fast modes which are usually decaying modes, 
allowing the simple strategy of removing the modes when they have 
become negligible. In the general case, one may encounter fast modes 
which are highly oscillatory and lightly damped. For such “WI\B” 
type problems, the current CSP algorithm does not work. In fact, for 
this type of non-linear problems, it is most likely that the system will 
behave stochastically in the slow time scale. 

During the course of the present research, the so-called RNG theory of 
turbulence made its appearance in the literature. In the RNG theory, 
the highly oscillatory and lightly damped modes of turbulence are for- 
mally “removed” from the solution, and the collective effects of these 
fast modes are emulated by an eddy viscosity. Using concepts and in- 
sights gained from the CSP theory, a critique of the RNG theory 1 
was published in 1992. This paper clarified some of the controversial 
procedures used in the RNG literature, and provided an alternative 
viewpoint. Most importantly, it provided an explanation for the out- 
standing quantitative successes of the original RNG theory. 

The Grant supported graduate students Mr. David Konopka and Ms. Xin 
Zhu who did their master degree theses on CSP, and Mr. Andrew Tron who 
is currently doing his Ph.D. thesis on CSP. In addition, the Grant supported 



Final Report for NASA Grant NAG-1-726 


12 


Dr. G. A. Goussis as a research staff member who now continues his CSP 
research as a faculty member of the University of Patra in Greece. 

Research on CSP will continue at Princeton under a AFOSR URI grant 
on aerothermochemistry which commences in 1993. 

6 Publications and Presentations 

Refereed Publications 

S. H. Lam and D.A. Goussis; “Basic Theory and Demonstration of Com- 
putational Singular Perturbation for Stiff Equations, Modelling and Simu- 
lation of Systems, P. Breedveld et. al. Eds., J.C. Baltzer A. G. Scientific 
Publishing Co., pp. 303-306, 1989. 

S. H. Lam and D.A. Goussis; “Understanding Complex Chemical Kinetics 
with Computational Singular Perturbation,” The 22nd Symposium (Interna- 
tional) on Combustion, pp. 931-941, August 1988. 

S. H. Lam and D. A. Goussis; “Conventional Asymptotics and Computa- 
tional Singular Perturbation for Simplified Kinetics Modelling,” Chapter 10 
of Lecture Notes in Physics, Reduced Kinetic Mechanisms and Asymptotics 
Aproximations for Methane Air Flames, M. Smooke, Ed., Springer- Verlag, 
1991. 

S. H. Lam and D. A. Goussis; “Sensitivity Analysis of Complex Simula- 
tions Using Basis Vectors,” Proceedings of the 13th IMACS World Congress 
on Computation and Applied Mathematics , R. Vichnevetsky and J. J. H. 
Miller, Eds, Volume 3, p. 1992, Criterion Press, Dublin, Irland, 1991. 

S. H. Lam; “On the RNG Theory of Turbulence,” Physics of Fluids A, 4 , 
5, May 1992. 

D. A. Goussis and S.H. Lam; “A Study of Homogeneous Methanol Oxida- 
tion Kinetics Using CSP,” 24th Symposium (International) on Combustion, 
The Combustion Institute, Pittsburgh, PA, pp. 113-120, 1992. See Appendix 
II. 

S. H. Lam; “Using CSP to Understand Complex Chemical Kinetics,” 
Combustion Science and Technology, 89, pp. 375-404, 1993. See Appendix 
I. 
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Conference Presentations 

S. H. Lam, D. A. Goussis and D. Konopka; “Computational Singular Per- 
turbation for Reacting Gasdynamics, Numerical Demonstrations I,” Open 
Forum Session, AIAA 8th Computational Fluid Dynamics Conference, Hon- 
olulu, Hawaii, June, 1987. 

S. H. Lam and D. A. Goussis; “Derivation of Simplified Multi-Step Reac- 
tion Models Using Computational Singular Perturbation,” Combined Meet- 
ing of The Combustion Institute (Eastern Section) and the National Bureau 
of Standard’s Annual Conference on Fire Research, Gaithersburg, Maryland, 
November 2-6, 1987. 

S. H. Lam and X. Zhu; “The Split CSP Method for Reacting Flow Prob- 
lems,” presented at the 41st Annual Meeting of the Division of Fluid Dy- 
namics of the American Physical Society, SUNY Buffalo, NY, November, 

1988. 

S. H. Lam and D. A. Goussis; “Construction of Reduced Chemical Kinet- 
ic Mechanism for a Methane Air Reaction System Using the CSP Method,” 
presented at the Reduced Kinetic Mechanism and Asymptotic Approxima- 
tions for Methane-Air Flames Workshop at UCSD, La Jolla, Ca., March, 

1989. 

S. H. Lam, D. A. Goussis and D. Konopka; “Time- Resolved Simplified 
Chemical Kinetics Modelling Using Computational Singular Perturbation,” 
AIAA 89-0575 Aerospace Sciences Meeting, Reno, Nevada, January 9-12, 
1989. 

S. H. Lam and D. A. Goussis; “Reduced Chemical Kinetics for Hydro- 
gen - Oxygen Reaction Systems,” Proceedings of the 19th Spring Technical 
Meeting, The Combustion Institute, pp. 121-127, Dearborn Michigan, April 
30-May 2, 1989. 

S. H. Lam and D. A. Goussis; “On the Batched Catalytic Fluidized Beds 
Problem,” presented at the 42nd Annual Meeting of the Division of Fluid 
Dynamics of the American Physical Society, NASA Ames Research Center, 
November, 1989. 

S. H. Lam, D. A. Goussis and P. A. Gnoffo; “Reduced and Simplified 
Chemical Kinetics for Air Dissociation Using Computational Singular Per- 
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turbation,” presented at the AIAA 28th Aerospace Sciences Meeting AIAA 
90-0644, Reno, Nevada, January 9-12, 1990. 

S. H. Lam and D. A. Goussis; “Computational Derivation of Simplified 
Kinetics Models,” presented at the 23rd International Conference on Com- 
bustion, University of Orleans, France, July 22-27, 1990. 

S. H. Lam and D. A. Goussis; “Sensitivity Analysis of Complex Chemical 
Systems,” presented at the 1990 Reduced Kinetic Mechanism and Asymptot- 
ic Approximations for Methane-Air Flames Workshop, Cambridge, England, 
July, 1990. 

S. H. Lam and D. A. Goussis; “Sensitivity Analysis of Complex Simula- 
tions Using Basis Vectors,” presented at the 13th IMACS World Congress 
on Computation and Applied Mathematics, Dublin, Irland, July 22-26, 1991. 
Also available as Report #1897-MAE, Princeton University, July, 1990. 

S. H. Lam; “The Effects of Fast Chemical Reactions on Mass Diffusion,” 
presented at the Annual APS Division of Fluid Dynamics Meeting at Talla- 
hassee, Florida, November, 1992. 

Unpublished Reports 

S. H. Lam; “Computational Singular Perturbation Procedure for Stiff 
Equations,” Mechanical and Aerospace Engineering Report No. 1697-MAE, 
Princeton University, 1985. 

S. H. Lam; “On Steady-State and Partial Equilibrium Approximations 
for Chain Reactions,” Mechanical and Aerospace Engineering Report No. 
1722-MAE, Princeton University, 1985. 

S. H. Lam and D. A. Goussis; “Research on Computational Singular Per- 
turbation (Progress Report $1: Basic Theory,” Mechanical and Aerospace 
Engineering Report No. 1779(a)-MAE, Princeton University, 1987. 

S. H. Lam and D. A. Goussis; “Research on Computational Singular Per- 
turbation (Progress Report $2: Numerical Demonstrations), Mechanical 
and Aerospace Engineering Report No. 1779(b)-MAE, Princeton University, 
1987. 

S. H. Lam and D. A. Goussis; “Understanding Complex Chemical Kinet- 
ics with Computational Singular Perturbation, Mechanical and Aerospace 
Engineering Report No. 1799-MAE, Princeton University, 1988. 
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S. H. Lam and D. A. Goussis; “Recent Developments in Simplified Ki- 
netics Modelling for Large Reaction Systems,” Princeton University Report 
#1864-MAE, 1989. 

Work in Progress 

S. H. Lam; “On Reacting Flows with Diffusion and Fast Chemical Reac- 
tions,” MAE Report #1999, 1992. See Appendix III. 

S. H. Lam and D. A. Goussis; “The CSP Method for Simplifying Kinet- 
ics,” MAE Report #1946, 1992. See Appendix IV. 

7 Epilogue 

With the support of this NASA Grant, we have been able to establish a firm 
mathematical foundation for CSP. The “correctness” of the results is beyond 
dispute — the numerical solutions of the full-model and the CSP-derived sim- 
plified model are guaranteed to stay below a user-specified threshold theo- 
retically, and confirmed by the computer code CSP91. The idea that com- 
putational power can be exploited to yield physical insights is new and will 
surely become more important in the future. 
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Abstract 

The conventional methods of simplified kinetics modeling through 
the use of partial- equilibrium and quasi-steady approximations are re- 
viewed and critiqued. The method of computational singular pcrtuv~ 
bation (CSP) is then presented with special emphasis on the interpre- 
tation of CSP data to obtain physical insights on massively complex 
reaction systems. A simple example is used to demonstrate how CSP 
deals with complex chemical kinetics problems without the benefits of 
intuition and experience. 


1 Introduction 

An ideal scenario in the (near) future for the study of chemical kinetics 
would be that a comprehensive, reliable and up-to-date database of vali- 
dated reaction rates is readily available to any researcher interested in any 
reasonable reaction system of interest. Using a suitable stiff solver, one can 
routinely compute for the numerical solution of a massively complex chem- 
ical kinetics system. However, the extraction of physical insights about the 

•This work is supported by NASA Langley’s Aerothermodynamics Branch, Space Sys- 
tems Division. 
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reaction system from the massive printouts is a different matter, and remains 
a formidable task. Generally speaking, answers to the following questions are 
highly valued: 

1. How can a reduced reaction system , involving a much smaller set of 
chemical species and elementary reactions, be derived which can ap- 
proximate the original massive reaction system both qualitatively and 
quantitatively? 

2. How can the rate-controlling elementary reactions be identified? What 
rate constants must be known accurately? What rate constants need 
not be known accurately provided that they are “sufficiently large”? 

3. Which elementary reactions can be considered fast reactions in the time 
interval of interest so that appropriate approximations can be applied 
to obtain useful approximate algebraic relations between the species? 

4. Which chemical species can be considered as “intermediaries” or “rad- 
icals?” in the time interval of interest so that their concentrations can 
be obtained from these approximate algebraic relations? 

In the present paper, we shall call approximate algebraic relations be- 
tween the species “ equations of state.” In addition, we shall use the word 
“radical” to connote a special meaning in the CSP context instead of its 
normal chemical structure context. An operational definition for a CSP rad- 
ical based on its mathematical role in the CSP theory will be given later. 
In most cases, a CSP radical will be found to be a chemical radical, and 
vice versa , but not always. To avoid confusion, we shall not use the term 
“intermediaries” altogether. 

Item $1 is important pragmatically because computational cost is a 
strong function of the number of chemical species. Questions in item # 2 
needs no justification; one can claim insights on the reaction system only if 
these questions can be satisfactorily answered. The conventional methodolo- 
gy deals with items #3 and #4 first, mainly by guessing based on experience 
and intuition. Some talented investigators simply know which reactions are 
fast and which chemical species are radicals under what conditions for certain 
reaction systems (Peters, 1985, 1991; Peters and Williams, 1987; Chelliah 
and Williams, 1990; Bilger et. al. 1990, 1991). For the less gifted, data from 
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numerical solutions of the full kinetics equations can be used to provide hints 
in the identifications. Once the fast reactions and the radicals are correctly 
identified, standard techniques are available to obtain answers to questions 
in items #1 and #2. 

The method of computational singular perturbation (CSP) (Lam, 1985; 
Lam et. al., 1988a, 1988b, 1989, 1991a, 1991b; Goussis et. al., 1990, 1992) 
provides a progarmmable mathematical algorithm to proceed routinely with- 
out the benefits of experience and intuition. It can be used to verified the 
validity of simplified models derived by ad hoc methods, and it can be used 
to deal with massively complex problems beyond the reach of such methods. 
Physical understanding of the reaction system under investigation can easily 
be extracted from the CSP data generated. 


2 Statement of the Problem 

Consider a reaction system of N unknowns 1 denoted by the column state 
vector y = [y 1 , y 2 , . . . , y N } T . The governing system of ODE is: 

I - «<*> (1) 

where 

g(y) = £s r F r (y), ( 2 ) 

r=l 

R is the number of elementary reactions being included in the reaction sys- 
tem, s r and F r (y) are the stoichiometric vector and the reaction rate of the 
r-th elementary reactions, respectively. The jV-dimensional column vector g 
is the overall reaction rate vector, and can be interpreted as the “velocity 
vector” of y in the iV- dimensional y- space. For a massively complex prob- 
lem, N and R can be large numbers, and the accuracy or reliability of the 
available rate constants is usually less than ideal. Assuming that adequate 
computing power is available, the computation of the numerical solution of 
(1) is not an issue. The challenge is to obtain answers to questions in items 
#1 and #2 in §1. 

iThe ^-dimensional column vector may include temperature, total density, etc. in 
addition to chemical species 
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As previously indicated, the conventional method relies completely on 
the skill of the investigator to identify the fast reactions and the radicals. 
Once the identifications are made, the subsequent theoretical development 
proceeds relatively routinely. The interested reader is referred to Appendix B 
of Williams (1985) for a detailed exposition of the conventional methodology. 


2.1 The Conventional Methodology by Example 

We shall use a simple hypothetical reaction system to demonstrate the con- 
ventional methodology. Let the state vector be y — [A, B , C\ where A and 
B are chemical concentrations and C is temperature. The reaction system 
consists of two elementary reactions: 

reaction# 1 : A + A *=* B, (3a) 

reaction# 2 : A ▼=* B. (3b) 


The stoichiometric vectors and the reaction rates are: 

Sl = [-2, 1, A H 1 ] t , F 1 = k x (A 2 — KiB), (4a) 

s 2 = [-1, 1, A H 2 } T , F 2 = k 2 (A — K 2 B). (4b) 


where the reaction rate coefficients k x , k 2 and the equilibrium constants K \ , I< 2 
are known and for the sake of simplicity their dependence on C is assumed 
negligible. The heat of reaction (in the proper unit) for the two reactions are 
denoted by AH X and A H 2 , respectively. 

The system of ODE is: 

-jr = SiF 1 + s 2 F 2 (5) 

dz 


which can be written out as follows: 


dA 

dt 

dB 

dt 

dC_ 

dt 


= - 2F 1 -F 2 , 

= F 1 + F 2 , 

= AH X F' + AH 2 F 2 . 


(6a) 

(6b) 

(6c) 
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The initial conditions, ,4(0), H(0) and C(0), are given. It can easily be 
shown that as t -* oo, A -> K x /K 2 and B -» K x jK\. Normalizing the 
unknowns intelligently, we have y» = [A», B., C.] T : 

A. = AKt/Kt, (7a) 

B. = BKl/I<u (7b) 

C. = CK 2 /(AH 2 K X ). (7c) 


The system of ODE for A*, B» and C, is: 

% = -4-KMAl - B.) - h(A. - B.), (8a) 

at A 2 

^ = I<MAl-B.) + K 2 k 2 (A.-B.), (8b) 

at 

+ (8c) 

at A 2A/22 

The initial conditions are: 

>1.(0) = A{0)K 2 /K u (9a) 

£.(0) = B(0)Kl/K u (9b) 

C.(0) = C(0)K 2 /(AH 2 K x ). (9c) 


As t —* 00 , we have A, — ► 1 and B, — * 1 . 

Because this reaction system is a hypothetical one, the concept of con- 
servation of atomic species cannot be applied — we do not know what atomic 
species are involved. However, it can easily be verified mathematically that: 

(AH X - A H 2 )A A (AH x - 2A H 2 )B + C = Constant (10) 

is an exact “integral of motion” valid for all time, and can be physically 
interpreted as a statement of conservation of total energy. 

For this simple example, if the “speeds” of the two reactions are competi- 
tive, there is no simplification available and little or no general statement can 
be made. However, when one reaction is much faster — in some sense — than 
the other, then certain mathematical simplifications are available and certain 
useful general statements can be made. Intuitively, we expect a rapid tran- 
sient period in which the fast reaction would dominate , followed by a slower 
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evolutionary period in which the slow reaction would be controlling. In gen- 
eral, we expect to be able to neglect some of the slower reactions in the rapid 
transient period. However, it is well known that the handling of the nearly 
exhausted fast reactions in the slow evolutionary period requires considerable 
care _it suffices to say here that they cannot be summarily neglected. 


2.2 The Case of reaction #1 being faster << 1) 

In the rapid transient period, the fast reaction #1 dominates. Hence, the 
slow reaction #2 can be neglected to yield the simplified model: 


dA, 

dt 

dB. 

dt 

dC, 

dt 


A 2 

I<MAl ~ &•), 


k 2 ah 2 


kmaI-b.). 


(lla) 

(lib) 

(11c) 


with (9a), (9b) and (9c) serving as the initial conditions. It is easy to show 
that in addition to (10), the following 


K 2 A « ■+■ 25, ~ 


^(A(0) + 2B(0)) 

A i 


( 12 ) 


is approximately “conserved” (remains constant) during this period. 

At the end of this period, the fast reaction #1 will have become “ex- 
hausted,” and A . and B , will be approximately related by B , « A\. In the 
slow evolutionary period which follows, an approximate equation of state is 
obtained by setting F 1 to zero. This is called the partial- equilibrium approx- 
imation and a special follow-up procedure is recommended. First, one of 
the participants of reaction #1 is somehow declared a radical. If F. is so 
declared, (8b) is used to eliminate the contribution from reaction #1 (i.e. 
K x k x {Al - F„)) from the rest of the equations. Then, the approximate e- 
quation of state is differentiated with respect to time and used to eliminate 
dB*/dt. The following simplified model is obtained: 


F. 

dA. 

dt 


Al 


K 2 k 2 

K 2 + 4 A. 


(A. - B.), 


(13a) 

(13b) 
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d C. 

dt 


[ 1 - 


Afh K 2 + 2A. 
AHffK 2 + 4A. 


)]h(A . - B .). 


(13c) 


It can easily be verified that the same answer is obtained (in contrast to 
results in §2.3.2 later) if A . is chosen as a radical. It is important to note 
that (13b) and (13c) are not obtained by substituting (13a) into (8a) and 
(8c). In fact, (13a) must never be used in (8a), (8b) or (8c). The initial 
conditions can be found from B* & A\ and (12): 


. , n ^ k 2 r sim+mm ,, 

M o+)« T [yi + 1J 


(14) 


with i?»(0+) and C,(0+) given by (13a) and (10). Note that ^ does not 
appear at all in (13a), (13b) or (13c), hence its value must be unimportant 
provided that it is “sufficiently large.” In this time period, (13a), (13b) and 
(13c) together behave as a “one-step” reaction with an effective stoichiomet- 
ric vector and an effective reaction rate , and generate approximate solutions 
which satisfy (10) exactly. 


2.3 The Case of reaction #2 being faster << 1) 

In the rapid transient period, reaction #1 is neglected to yield the simplified 
model: 

f» —k 2 (A* — B m ), (15a) 

« K 2 k 2 (A* - B m ), (15b) 

« k 2 (A m — B r ). (15c) 

again with (9a), (9b) and (9c) serving as the initial conditions. It is easy to 
show that 2 

K,A. + B. « fy(A( 0) + B( 0)) (16) 

is approximately conserved (remains constant) during this period. 

At the end of this period, A„ and B, will be approximately related by 
A*. For this case, we shall present both the the partial-equilibrium 
approximation and the quasi-steady approximation in the slow evolutionary 
period. We shall see that they produce different answers. 


dA. 

dt 

dB, 

dt 

dC. 

dt 
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2.3.1 The Partial-Equilibrium Approximation 

In the slow evolutionary period, we can proceed as before by assuming el- 
ementary reaction $2 to be in partial-equilibrium. Using the same recom- 
mended procedure, we obtain the following simplified model: 


B. « A„ 

(17a) 

dA . K\k\ 2n 

dt K 1 + k 2 {b - A - } ’ 

(17b) 

dC . ( I<2 + 2 AH x ,I< x k x , D ^ 

dt a l K, + 1 A H, ] K, {B - ' 4J - 

(17c) 

with initial conditions: 


A(0+)»* MO) + B(0)) 
A.(U+) ~ Kx (1 + K 2 ) 

(18) 


with B.{ 0+) and (7.(0+) given by (17a) and (10). It is emphasized again that 
(17b) and (17c) are not obtained by substituting (17a) into (8b) and (8c). 
It can be shown that the partial-equilibrium approximation above can be 
formally justified when K x k x << K 2 k 2 and A»(0) = 0(1), with no restriction 
on the magnitude of K 2 , and the resulting approximate solutions satisfy (10) 
exactly. 

2.3.2 The Quasi-Steady Approximation 

Alternatively, we can declare either the species B * or A, to be a radical 
in this period. The quasi-steady approximation neglects the time derivative 
term from the ODE of the radical to yield an approximate equation of state. 
If B . is so declared, the approximate equation of state obtained from (8b) 
is used to simplify (8a) and (8c). The resulting set of equations (without 
(8b)) is then augmented by the approximate equation of state to yield the 
simplified model. We obtain: 


B. 


(19a) 

dA. 

dt 

K\k\ ( j 2 \ 

~ jr 

(19b) 

dC. 

dt 

n AH x ^K x k x 2 \ 

* 11 - A H 2 ] K , (B - 

(19c) 
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with initial condition (18). Note that (19b) and (19c) agree with (17b) and 
(17c) and satisfy (10) only when K 2 >> 1. If K 2 « 1, then A , should have 
been declared a radical instead, and a different result valid only in this limit 
is obtained. If K 2 is of 0(1), then the quasi-steady approximation is simply 
incorrect — while the partial-equilibrium approximation is valid for all K 2 . 

Again, the value of k 2 does not appear in the simplified model in the slow 
evolutionary period. The inference is that its value is unimportant provided 
that it is sufficiently large. 

2.4 Comments on Conventional Asymptotics 

The above example prompts the following questions and observations: 

1. Exact algebraic “conservation laws” such as conservation of atomic 
species and conservation of energy frequently exist in chemical kinetics 
problems. In the above example, (10) is such a conservation law. Can 
such conservation laws, when they exist, be mathematically identified 
and derived? 

2. In the rapid transient period, additional temporary approximate con- 
servation laws may exist. In the above example, (12) and (16) are 
temporary approximate conservation laws valid for the respective cas- 
es. Can such temporary approximate conservation laws, when they 
exist, also be mathematically derived? 

3. When a fast reaction exhausts itself, its net reaction rate is not neg- 
ligible in general from the governing equations. For example, when 
reaction #1 is the faster reaction, the exhausted F 1 is not zero in the 
slow evolutionary period but is, to “leading order,” given by: 



KihiAl-B.) - 

(K 2 + 2A„) . g , 

(A 2 + 4A.) 2 2( * 

(20a) 

or 

F 1 

K X +2A 2 
tfi+4A ’ 

(20b) 


which can be derived by comparing (8a) with (13b). In other words, 
when a fast elementary reaction is exhausted, its net reaction rate may 
be competitive with the slower reactions. 
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4. In any slow evolutionary period, the temporary approximate conser- 
vation laws valid in the previous periods are replaced by approximate 
equations of state obtained by applying the partial-equilibrium or the 
quasi-steady approximation. In our example, (13a), (17a) and (19a) are 
such equations of state for the respective cases. Can these equations of 
state be mathematically derived? How is (13a) related to (20b)? 

5. In a slow evolutionary period, the values of the exhausted fast reaction 
rate constants are not important provided that they are fast enough. 
In our example, the fastest rate coefficient (either or k 2 ) does not 
appear in the simplified models when M = 1. From the computation 
point of view, their absence means the simplified model is no longer 
“stiff” (Aiken, 1985). 

6. To get started on the conventional method, one must somehow decide 
that certain reactions are faster than others based on some vague, in- 
tuitive judgment. Unlike linear problems, the “speed” of a “reaction 
mode” in a non-linear problem can depend strongly on the state of the 
system. For example, even if K 2 k 2 >> K\ki, reaction #1 can still 
dominate reaction #2 initially provided that A,( 0) is sufficiently large. 
Since we are clearly interested in non-linear problems, a mathematical 
criterion for the identification of fast reactions is very much needed. 

7. Even after the fast reactions have somehow been identified, there is 
still the question on whether the partial-equilibrium or the quasi-steady 
approximation should be used to generate the approximate equations 
of state, and which species — radicals — can be eliminated by using these 
algebraic equations. 

8. In Williams’ book (1985), a radical is described as a chemical species 
which is “neither initial reactants nor principal products,” a description 
too vague to be useful. In §2.3, we saw that whether A m or B m could 
justifiably be declared a radical depends on the magnitude of K 2 , and 
that when K 2 = 0(1) the quasi-steady approximation is simply wrong. 

9. We shall adopt a pragmatic but precise definition for a radical in the 
CSP context: a chemical species is a CSP-radical if its ODE can be 
replaced by an approximate equation of state. For the moment, we 
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shall leave the questions of how to identify the radical(s) and how to 
solve for them from CSP-derived equations of state until the concepts 
of “radical pointer” and “radical correction” are introduced later. 

10. We emphasize here again that the approximate equations of state ob- 
tained from the partial-equilibrium approximation must never be used 
directly in the original governing equations — the so-called recommend- 
ed procedure must be followed. However, the approximate equations 
of state obtained from the quasi-steady approximation may be used 
more freely, but the validity of the results is more restricted. In fact, 
the quasi-steady approximation can be viewed as an additional ad hoc 
restriction imposed on the results of an appropriately applied partial- 
equilibrium approximation. 

11. Conventional asymptotics assumes that the investigator is not only ca- 
pable of non-dimensionalizing the variables intelligently, but also com- 
petent to take advantage of the presence of large or small parameter(s) 
in the resulting formulation. Can a theoretical structure be developed 
such that intelligent non-dimensionalization is not required? 

12. In general, the conventional “magnitude” of a dimensional vector y 
computed using the standard inner product may not make physical 
sense. In the above example, ||y|| 2 = ( A 2 + B 2 + C 2 ) 1 / 2 makes no 
physical sense at all. In the language of mathematics, the “norm” of 
the vector space of y needs to be defined. Hence if intelligent non- 
dimensionalization is not assumed, the concept of a “good approxima- 
tion” needs to be explicitly clarified. 

3 The Goal and the Idea of CSP 

The conventional method is only viable for relatively simple problems for 
which adequate amount of experience and intuition have been accumulat- 
ed, and that the algebra involved is manageable. For massively complex 
problems with large values of N and R, a better method is clearly needed. 
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3.1 The Goal 

The goal of CSP is to computationally derive time-resolved simplified models 
for massively complex chemical systems assuming that the appropriate reac- 
tion rate database is available. The novel feature of CSP is that the power 
of the modern computer is exploited not only to construct the numerical 
solutions, but also to derive the approximate equations. 

The output of CSP is a set of numbers called the CSP data. An essential 
attribute of CSP data is that it must be easy to interpret. Answers to 
questions such as: 

• how many approximate equations of state are available and how can 
they be derived? 

• which species can be considered radicals — in the CSP context so that 
they can be computed from these equations of state? 

• which elementary reactions are controlling the reaction rate of the sys- 
tem? and 

• what is the minimum reduced reaction system that will generate an 
approximate solution with a user-specified accuracy? 

and others will be provided by the CSP data. 


3.2 The Idea 

The physical problem is completely specified by g(y), a non-linear function 
of y obtained by summing all the physical processes which contribute to the 
time rate of change of y. The question is: is there a better representation 
than the physical representation? 

A representation is said to be a physical representation when the theoreti- 
cian formulating the problem can explain and interpret each term physically. 
Equation (2) is such a representation because as written each term represents 
the contribution of each of the R reactions included in the full model. Must 
we always write g in this particular form? Since g is a iV-dimensional vector, 
it can always be expressed in terms of a set of arbitrarily chosen N linearly 
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independent column basis vectors, a it i = l,2,---,N. When some other set 
of basis vectors is used, g has an alternative representation : 

g = £>,/ (21) 

*=1 

where /*, called the amplitude of the i-th mode, is given by: 

/'(y)sb'Gg = f;B;F', . = 1,2,. ..,1V, (22) 


where 


B' r = b‘ © s r , i = 1,2, ... , JV, r = 1,2,...,#, 


(23) 


and © is the dot product operator of the iV-dimensional vector space. When 
properly normalized, f' can also be considered as a “progress variable” for 
the i-th mode. The set of N row vectors b* are the inverses of a,; together 
they satisfy the following orthonormal condition: 


b* © a_, = 8), i,j 


1,2,. ..,N. 


(24) 


Note that once the set a, is chosen, the associated set b‘ is straightforwardly 
computed. Note also that B' r is not necessarily dimensionless unless deliber- 
ately made so. 

In our example, we have N = 3. The physical representation chooses (by 
default) the following column basis vectors: 


a x = s x , 


(25a) 

a 2 = s 2 , 


(25b) 

a 3 = s 3 = [0, 0 

1] T . 

(25c) 


where we have added S3 to form a complete set. Any vector which is linearly 
independent of a x and a 2 may also be used. The associated row basis vectors 
are obtained by solving (24): 


b 1 = [-1, -1, 0], 

b 2 = [1, 2, 0], 

b 3 = [Atfi-Atf 2 , AH 1 -2AH i , 1]. 


(26a) 

(26b) 

(26c) 
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We can readily verify that: 

Z 1 = b 1 ©g = J F 1 . ( 27a ) 

/ 2 = b 2 © g = F 2 . (27b) 

f 3 = b 3 0 g = F 3 = 0. (27c) 

The CSP idea is very simple: instead of using the physically meaningful 
stoichiometric vectors as the default basis vectors, let’s exploit the theoreti- 
cians’ prerogative of trying different alternatives — may be something else 
works better. 


3.3 Same Problem, Different Basis Vectors 

We shall rework the first case of our example using dimensional formulation. 

Instead of the physically meaningful stoichiometric vectors, we shall use, 
without comments at this point, the following set of basis vectors: 


a x = [-2, AH 2 ] 


I T _ 


Si, 


(28a) 


a 2 = {K 1 ,2A,(Ki+4A)AH 2 -(Ki+ 2A)AH 1 ] t , (28b) 

2 K 1 +4A [ ’ v 

a 3 = [0, 0, If. ( 28c ) 

The following row vectors are their inverses and together they satisfy (24): 

1 


b ‘ = k^Ta [ - 2A ' k " 0| ' 

b 2 = [1, 2, 0], 

b 3 = [AHi — AH 2 , AHi — 2AH 2 , 1]. 


7C x +2A f2 

K\ -j- 4A 


The amplitudes of the modes are: 

f 1 = b 1 ©g = F 1 + 

/ 2 = b 2 © g = F 2 , 
f 3 = b 3 O g = F 3 = 0. 

In terms of these basis vectors, the original reaction system becomes. 

dv 


(29a) 

(29b) 

(29c) 

(30a) 

(30b) 

(30c) 




i 


•2 


•3 
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We can rewrite (31) in long-hand notation as follows: 


dA 

dt 

dB_ 

dt 

dC_ 

dt 


-2 f 1 + 
f' + 


Ki 

Ki + 4A 

2A 




Ki+AA 


f\ 


Atf,]/' + [A^ - 


(32a) 

(32b) 

(32c) 


Taking the inner product of b’ with (31), we obtain: 


1 


d 


!<! + 4 A dt 


—[(AHi - A H 2 )A + (A H t - 2A H 2 )B + C] 
dt 


[-A' + K.B] = f\ 
\aa 2B) = f 2 , 
f 3 - 


dt 


(33a) 

(33b) 

(33c) 


We shall assume for the moment that somehow the modes have been ordered 
in descending speed — mode $1 is faster than mode $2. This is physically 
equivalent to saying that F 1 is estimated to be faster than F 2 . Reaction rate 
/ 3 of mode #3 is, for this problem, identically zero. 


3.3.1 Exact Conservation Law 

Whenever the amplitude of a slow mode is considered “negligible,” it can be 
summarily neglected. Since for this problem / 3 = 0, mode #3 is always neg- 
ligible, (33c) recovers (10), the conservation law for energy obtained earlier. 
However, a dormant mode does not always yield a conservation law. See §7.2 
later. 


3.3.2 The Rapid Transient Period 

In the rapid transient period, we assume that the dominant mode is f , while 
f 2 and / 3 are negligible in comparison. 3 The simplified model is: 

< 34 > 

2 Compare them with (6a), (6b) and (6c). 

3 A slower mode may not necessarily be negligible in this period; if it is not, it must be 
kept. 
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Neglecting the right-hand side of (33b) and (33c), we recover (12) and (10), 
the temporary and permanent conservation laws obtained earlier. The so- 
lution of (34) automatically satisfies these constraints. As time progresses, 
mode #1 decays rapidly and becomes exhausted as p — ► 0, and y approaches 
a fixed point. 

3.3.3 The Slow Evolutionary Period 

In the slow evolutionary period, (12) ceases to be valid, but the exhausted 
mode ^1 replaces it with a new approximate equation of state. 

/'( y) = EB,'F r (y)«0, (35a) 

r= 1 

= + < 35b) 
I<1 + 4A 

which is the counterpart of (13a), and agrees with it in the |F J | >> |F 2 | 
limit. However, unlike (13a), which we emphasized must never be substituted 
directly into the original governing equations, (35a) needs no such admonition 
and can be so used freely. Using it in (31), we obtain the simplified model 
for this period: 

% K 82/2 • (36) 

Solution of (36) automatically satisfies (35a) if it is satisfied initially. No 
decision on which species should be eliminated (t.e. considered as a radical) 
was needed so far. All results obtained can be shown to be consistent with 
the conventional method in the asymptotic limit of reaction #1 being much 
faster than reaction #2. 

3.3.4 Response to a Slow Third Reaction 

For our simple example with only two reactions and three unknowns, the 
reaction system reaches its steady state after both modes $1 and $2 are 
exhausted. What happens if there is a third reaction? Denote the stoichio- 
metric vector and reaction rate of the third reaction by S3 and F 3 respectively. 
In addition, assume that the three stoichiometric vectors are linearly inde- 
pendent, and that |F J | » |F 2 | » |F 3 | so that the added reaction is the 
slowest. 
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Using the same basis vectors in §3.3, we can proceed as before. The main 
impact of this complication is that /’ = b’ Og now contain contribution from 
F 3 . After mode #2 is exhausted, a even slower mode #3 will take over. One 
can easily verify that the simplified model for this period is: 

% * as/3 - (37) 

Note that in this time period, the response of A and B to the slow third 
reaction is negligible; only C responds significantly. 

3.4 The Obvious Questions 

It can be shown that the three simplified models (34), (36) and (37) dis- 
played above are the leading-order approximations in the asymptotic limit of 
|F X | >> |F 2 | >> |F 3 | for the indicated time periods, as obtained earlier in 
§2.2. Hence, the alternative representation of g using the special set of basis 
vectors ((28a), (28b) and (28c)) is clearly a good idea in this limit. 

It can easily be verified that (28a), (28b) and (28c) are the right eigen- 
vectors of the Jacobian of g in the limit of k\K\ >> k 2 K 2 ■ Why do these 
special basis vectors allows us to summarily neglect the exhausted fast modes, 
while the same casualness would not be tolerated previously? What is the 
relationship between the conventional methods and the eigen-vectors of the 
Jacobian? How does one know that the basis vectors which worked well for 
one set of initial conditions will work for a different set of initial conditions? 
What happens when k 2 K 2 /kiKi is only moderately small? 

Additional questions are: How does one extract physically interesting 
information from these basis vectors? Which species can be considered a 
radical? How can the accuracy of these models be improved? How can the 
whole procedure be generalized to deal with a massively complex problem? 


4 The Ideal Basis Vectors 

Given any set of N linearly independent basis vectors, one can always de- 
compose the A^-dimensional vector g into N additive components or modes. 
Most investigators use the physically meaningful stoichiometric vectors as 
basis vectors by default. When the set used is non-ideal, the speed ranking 
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of the modes is unclear, and the residual net reaction rate of an exhausted 
fast reaction mode can not be neglected in general, and must be treated with 
great care. We have seen in the above example that the use of an intelligently 
chosen alternative set of basis vectors can make a difference. A good set of 
basis vectors not only should order the modes according to their speed, but 
also guarantee that the residual net reaction rate of an exhausted mode be 
negligible. 

To find out how the amplitude of each mode evolves with time, we differ- 
entiate (1) with respect to time to obtain: 

4=J0g (38) 

dt 

where 

J = — = N x N Jacobian matrix of g with respect to y. (39) 

dy 

Since the solution y(t) and therefore J(t) are both known functions of time, 
(38) can be considered a linear ODE (with known time-dependent coeffi- 
cients) for g. Taking the inner product of b‘ with (38) and using (21), we 
obtain: 

C = S>;/ J - i = 1.2 N, (40) 

at J=1 

where 

Aj = [~- + b' 0 J] O aj, i,j = 1,2, . . . , AT. (41) 

J dt 

A set of basis vectors a i{t) is said to be ideal if (i) the inverse row vectors 
b '(t) can be accurately computed from (24) for all time interval of interest, 
(ii) A'At) is diagonal, and (iii) the diagonal elements of A’ (t) are ordered in 
descending magnitudes. For linear problems where J is a constant matrix, 
the ideal basis vectors would be the (constant) ordered eigen-vectors of J. 
For non-linear problems, the eigen- vectors of J are time-dependent, and they 
do not diagonalize A* . 

5 The CSP Refinement of Basis Vectors 

The method of CSP does not attempt to find the ideal set of basis vectors 
even when g is linear. Instead, it assumes that, at any moment in time, a 
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trial set of ordered basis vectors is somehow available, that the first M fastest 
modes are exhausted as measured by some criterion (to be specified later), 
and generates from this trial set a new refined set of basis vectors, a? and 
b* i = 1 2 . . . ,N using a two-step refinement procedure (Lam and Goussis, 
1991a). When recursively applied, the refinement procedure successively 
weakens the coupling between the fast modes and the slow modes. 

The step #1 refinement is: 


b™(M) = b m + £ p7{M) b J , 

m = 1,2, ... , M, 

(42a) 

J=Af+ 1 



M 

a °j(M) = aj-^a n p3(M), 

J = M + l,...,iV. 

(42b) 


n=l 


The step #2 refinement is: 

M 

b'(M) s b' - Y I = M + 1,...,N, (43a) 

n=l 

a ° m (M) = a m + a m = l,2,...,M. (43b) 

J=A/+1 


The matrices pJ{M) and q J m {M) above are defined by: 


M 

p?(m) = £c(M)A3, 


n=l 

M 




m 


m 


1,2 

1,2,..., Af, 


and t”(M) is the inverse of A” (M): 

M M 

Y, r:(U)K k m (M) = y wl w = 

k = 1 k=l 


J = M + (44a) 

J = M + 1,...,JV, (44b) 

n,m = 1,2, ..., M. (45) 


The refinement procedure is readily programmable, and the refined basis 
vectors after each step satisfy the orthonormal relations, (24). It is essentially 
a generalization of the so-called “power method” for computing eigen- vectors, 
and produces a block diagonal A* when converged. 
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The amplitude of the m-th fast mode satisfies the following differential 
equation (Lam and Goussis, 1991a): 


where 


df m 

dt 


M 

E a”( r - iz). 


n=l 


m 


1,2 M, 


f: 

n 


= - e p”/ J , = 

J=M+1 


= b-©g = r-/z, 


1,2 ,...,M. 
m = 1,2, , M. 


(46a) 

(46b) 

(46c) 


Note that is some linear combination of the trial slow mode amplitudes. 
Solving for from (46c), we obtain: 


C = / m -/o m = (b m -bD©g, m = 1,2,..., M, 

which shows clearly that is proportional to the change of b m to b™ as a 
consequence of the step #1 refinement procedure. In general, /™, which is 
computed using the trial set of basis vectors, can not be considered “small.” 
Let r(M) denote the magnitude of the time scale of the slowest of the 
fast modes. Let t(M + 1) denote the magnitude of the currently active time 
scale. Theoretically, they can be estimated by the M-th and (M + l)-th 
eigen-values of J. In the limit of small r(M)/r(M + 1) — large time scale 
separation — the formal asymptotic solution of (46a) is: 

r = /» + /"<« + •■■. m = l,2 M, (47a) 

Co = m=l,2,...,M. (47b) 

n=l al 

The order of magnitude of can be estimated by: 

Co « 0(fZr(M)/T(M + 1)), m = 1, 2, . . . , M. (48) 

Assuming that r(M + 1) >> r(M), the order of magnitude of is now 
small in comparison to /” — i. e. it is a “higher order” term. We have 
adopted the notation that superscript or subscript o indicates a variable 
evaluated with refined basis vectors. The above derivation can be found in 
Lam and Goussis (1991a). 

Maas and Pope (1992) recognizes the importance of the eigen- vectors 
of J, but recommends the use of its Schur vectors (Noble, 1988) — which 
transform J into a lower- triangular form — as basis vectors. The advantage 
of Schur vectors for dealing with dimensional vectors is not clear. 
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5.1 The Trial Set of Basis Vectors 

To get started, CSP needs to have a trial set of basis vectors. A number of 
options are available. 

For sufficiently simple problems, one may wish to proceed analytically. 
N linearly independent stoichiometric vectors chosen from the reactions ac- 
tually included in the full model (with supplementary basis vectors added as 
required) may be used as the default trial basis vectors. Educated guesses 
are then needed to establish the speed ranking of the modes. 

A fool-proof procedure is to compute for the eigen-vectors of J at t = 0 
and use them as time-independent trial basis vectors for t >0. The recip- 
rocal of the eigen-values, denoted by r(i), is an approximate measure of the 
characteristic time scales of the modes. We shall assume that the t(z) s are 
essentially real and are ordered in ascending magnitudes. 

When the refinement process is performed numerically on a computer, the 
refined basis vectors used in the previous time-step can be used as the trial 
basis vectors for the new time-step. Under this strategy, the initial choice of 
trial basis vectors at t = 0 is not critical. 


5.2 The Number of Exhausted Modes 


The refinement procedure requires the knowledge of M , the number of ex- 
hausted modes. In general, none of the fast modes are negligible at t — 0 
(i.e. M = 0). As time progresses, the faster modes will eventually decay and 
become exhausted (t. e. f™ — * /™oo) - The number of exhausted fast modes, 
M(t), is determined by requiring that their contributions to g (see (54) and 
(55a) later) over the next time interval of 0(t(M + 1)) are negligible: 

K/ 0 "t(M)| < y error , m = 1,2, . . . M, (49a) 

K/oV^ + !)l < y error j m = 1, 2, . . . , M, (49b) 

where y erTO r is a dimensional column vector representing a user-specified 
threshold of absolute error allowed, and the vector inequality is individually 
enforced for every components of the vector. If one relaxes the requirement 
that the accuracy threshold for each component of y be individually enforced, 
then the following are acceptable alternatives to (49a) and (49b): 


1/71 < 0 ( 


b™ 0 Yerror 

r(M) 


) , m — 1,2,..., Af . 


(49c) 



Appendix I 


39 


LCJ < o( 


Koy e 


t(M + 1 ) 


), m = 1,2, , M. 


(49d) 


Either of the above exhaustion criterion can readily be implemented compu- 
tationally. The value of t(M + 1) can be estimated by the current integration 
step-size selected by any automatic (explicit) ODE integration package. Note 
that the value of M need not be a monotonic function of time. 

Actually, the exhaustion criterion for the so-called CSP radicals defined 
in §6.4 can be made more lenient when the so-called radical correction (to 
be discussed in §6.5 later) is applied. A full discussion of this subtle point, 
however, is beyond the scope of this paper. 

For most problems, the desired time resolution At is usually known. If 

At > t(M), (50) 


then the M fastest modes must be nearly exhausted in the At time scale 
of interest. If we are not interested in the details of what happens in the 
next t(M ) seconds, we may use the simplified model with a non-zero value 
for M starting at t « 0, provided that we also adjust the initial conditions 
appropriately as demonstrated in (14) and (18). See §6.5 later. 


5.3 The Refinement Process 

The CSP refinement procedure consists of the following two steps. 

1. Refinement of the fast row vectors using (42a) and the slow column 
vectors using (42b). This step depresses the magnitude of the upper- 
right block of A£, (m = 1, . . . , M, I< = M + 1, . . . , N ), by the factor 

+1), and thus weakens the coupling of the fast mode am- 
plitudes from the slow. 

2. Refinement of the slow row vectors using (43a) and the feist column 
vectors using (43b). This step depresses the magnitude of the lower- 
left block of A*, (K = M + 1, . . . , N, m = 1, . . . , M ), by the factor 

+ 1), and thus weakens the coupling of the slow mode am- 
plitudes from the fast. 

The full cycle of the two-step refinement process renders the new A’ 
matrix calculated from the refined basis vectors more nearly block diagonal 
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than before. Step #1 improves the accuracy of the simplified model, while 
step #2 ensures that the simplified model is not stiff. The speed ranking of 
the modes can be estimated by the diagonal elements of the A*- matrix. 

In our example in §3.3, the set of alternative basis vectors (28a), (28b) and 
(28c) used was obtained by adopting (25a), (25b) and (25c) as the ordered 
trial set, and then refining them analytically through step #1 using M = 1. 
We obtain: 


K = 

K t +4A 


1 

a 2 “ 

K,+4A 

a 3 = 

[ 0 , 0, 


[-2 A, Ki, 0] r , 

[K u 2 A, ( K 1 + 4A)AH 2 -(K 1 + 2A)AH l ] T , 

1] T 


(51a) 

(51b) 

(51c) 


These refined vectors are significantly modified from their trial counterparts. 
Performing step # 2 , we obtain: 




= (1 + 


_ 


2 Ae 


2 - 


I< it 


0 ], 


Ki + A A' " K\ + 4A' 

= \AH X -AH 2 , AH 1 -2AH 2 , 1], 

K\t . 2 At 


= [-2 + 


1 + 


K\ + 4A ’ " ' Kr + 4A' 

A//, - ( ^*^4 Aff, - A HM, 


(52a) 

(52b) 


(52c) 


where e is a dimensionless parameter defined by 


_ ^2(7^2 + 2 ) 
£= k)(Ki + 4A)' 


(53) 


These refined vectors are only slightly modified from their trial counterparts 
for c << 1. Generally speaking, whenever the set of trial a, basis vectors are 
intelligently chosen, step #2 will provide only small corrections. If the trial 
set was randomly chosen, more than one full cycle of the two-step refinement 
process may be necessary to generate the correct leading order result in the 
limit of t — * 0. 

The trial set is usually chosen initially to be time-independent for the 
sake of convenience. For non-linear problems, the refined set will in general 
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be time-dependent, and the associated time derivatives in Aj must be prop- 
erly evaluated. Analytically, this is a straightforward step which could be 
quite tedious algebraically. From the programming point of view, however, 
good approximation to db'/dt can be obtained in a number of ways, such as 
utilizing the availability of stored values in the previous time steps and/or 
the predicted values in the next time step used in most integration routines. 


6 Using the Refined Basis Vectors 


Using the refined basis vectors, the governing system of ODE’s become: 


^ = g 
dt * 


_ „o y jast 


+ S°o 


.o,slow 


(54) 


where 


M M 


,o t fast 

>0 

= E 

o tm 
m J o 

= ( E a m b 

r)©g, 

(55a) 


771 = 1 


771=1 



f m 

J o 

7 771 

= b 0 

Og, 

m = 1,2,.. 

.,M, 

(55b) 




R 



.OySlOW 
► 0 

= g - 

°,faat 

o o 

_ V' c°' slow 
2 —j S o,r 

r=l 

F, 

(55c) 



M 




.OySlOW 

'o y r 

= (I- 

■Ei 

;bD 0 s„ 

II 

to 

(55d) 


771 = 1 


The column vector is mathematically the projection of s r in the so- 

called slow subspace or manifold , and can be interpreted physically as the 
effective stoichiometric vector of the r-th elementary reaction. Unlike the 
original chemical s r , the components of s°f low are not necessarily integers or 
rational numbers, and may involve species which do not appear in F r ( y). 

When the M fast modes satisfy the exhaustion criterion decribed in the 
next section, the simplified model is simply: 


dy 

dt 


g 


OySlOW 

o 


N 

= E 


a° f K 
& KJo • 


K=M + 1 


(56a) 

(56b) 
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6.1 Classification of Modes 

In general, there are three kinds of modes: exhausted modes, currently active 
modes, and dormant modes. Let r(i) represent the characteristic time scale 
of the j-th mode. It is assumed that: 

T(m) << t(K), m = l,2,...,M, K = M + 1, . . . , N. (57) 

The ordering among T(m) and t(K) is not important. 

Exhausted Modes: By applying the exhaustion criterion (49b), the num- 
ber M of exhausted fast modes is determined, yielding a set of M 
approximate equations of state: 

f? = K 0 g = £ B Z r * 0, m = 1, 2, . . . , M, (58) 

r=l 


where 


B' ot = bj, 0 s r , i = 1,2,..., iV, r = 1,2,...,/?. 


(59) 


Currently Active Modes: The remaining N — M modes are kept in (56b). 
Using (55c), we rewrite the simplified model as follows: 



'tslfT 

r=l 


(60) 


Equation (60) can be numerically integrated without the need of a 
stiff solver. The integration time-step used should be a fraction of 
r(M + 1), the characteristic time scale of the (Af + l)th mode. The 
initial condition of (60) must satisfy (58) in the sense of (49c). If 
the initial trial fast column basis vectors were randomly chosen, the 
accuracy of (60) is + l)) c where c is the number of full 

two-step refinement cycles. If they were initially “intelligently” chosen, 
then c is the number of step #1 refinements. 

Dormant Modes: The amplitude of some of the slow modes kept in (56b) 
may be extremely small or even identically zero. When a slow mode is 
deemed negligible— when the its contribution to g°' alow in the current 
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time period is estimated to be less than y erroT component by compo- 
nent, it is said to be a dormant mode. In the example, mode #3 is a 
permanent dormant mode. If mode #N is a dormant mode, then: 

f? = 0 g » 0, (61) 

/„" = E K F ' ~ 0. (62) 

r=l 

However, unlike (58), (62) is not a useful equation — it can not be used 
as an approximate equation of state. Instead, (61) is the more useful 
equation, indicating that the vector g has negligible projection in a 
specific “direction.” If can be shown to be the gradient of a scalar — 
a condition known as holonomic constraint in classical mechanics — then 
that scalar is a conserved “integral of motion.” Conservation laws of 
atomic species manifests themselves as permanent dormant modes with 
constant row vectors (which are always holonomic). Equation (33c) is 
a conservation law since b 3 is holonomic and / 3 is identically zero. 

The solution of (60) automatically satisfy (58), (61) and (62). The 
distinction between (58) and (61) will be further discussed later. 


6.2 Equation of State and Participation Index 

The left-hand side of (58) usually consists of large positive and negative terms 
which nearly cancel each other. Which elementary reactions participate most 
strongly in this balancing act? The participation index, denoted by Pf a , is 
designed to provide this information and is defined as follows: 


pi = ^ 

T '° fl n.-H+ b * 0y “ 

l. p./ + znn 


, i = 1,2, ...,7V, r = 1,2, . . . , 


r(M +1) 


where it is assumed that forward and reverse reactions are counted as dis- 


tinct so that no cancellation occurs within any F r . Pf Q is a measure of the 
participation of the r-th elementary reaction to the balancing act of the i-th 
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mode. To get an idea of which elementary reactions are the major partic- 
ipants of each fast mode, we just need to sort the elementary reactions by 
their participation indices. A exhausted fast mode usually has several sig- 
nificant participation indices, a slow dormant mode may not have any, while 
a permanently dormant mode representing a conservation law has none. In 
our example worked out in §3.3, the participation index data for mode #1 
will show that reaction #1 is the main participant of / 1 as defined by (30a), 
reaction #2 is the main participant of / 2 as defined by (30b), while f 3 = 0 
hits no participant at all and can be interpreted as a conservation law. 


6.3 Simplified Model and Importance Index 

The right-hand side of (60) consists of terms which control the reaction rate 
of the system. For each component of the state vector y, we can order the 
terms on the right-hand side in descending magnitudes. Which elementary 
reactions are most important in controlling the reaction rate of a particular 
species of interest? The importance index, denoted by I' T '° and designed to 
provide this information, is defined as follows: 


7 - 1,0 

A o>r 


&i,o,slow J?T 


R 


s i,o, slower _|_ 


r=l 


*'= 1,2,..., AT, r = 1,2,...,#. 


r(M + 1) | 


(64) 


I' 0 ’° is a measure of the relative importance of the contribution of the r-th 
elementary reaction to the current reaction rate of the i-th species. Note that 
if /’’ s were computed using the default s r ’s, we would obtain the misleading 
information that the fast reactions are rate-controlling even after they are 
exhausted. In general, the values of the importance indices change discon- 
tinuous^ at the moment when a new fast mode is declared exhausted — the 
7*’®’s of the major participants in the new exhausted mode will drop, while 
those in the emerging rate-controlling modes will rise. In the example worked 
out in §3.3, the importance index data on (36) computed using s° 0 '* T low (see 
§7.5) will show that reactions #1 and #2 are important in controlling the 
overall reaction rate but each in their own time period only. 
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6.4 The Radical Pointer 


When M modes are exhausted, M algebraic equations are obtained from (58) 
which can be used to replace M of the ODE’s in the subsequent epochs. The 
question is: which ODE’s can be replaced? In other words, which species are 
the radicals? 

CSP associates with each exhausted mode one or more species by the 
radical pointer. The radical pointer of the m-th mode, Q m (i), is defined by 
the N diagonal elements of the N x N matrix a m b m , refined basis vectors 
preferred. Geometrically, the magnitude of Q m (i ), which is dimensionless and 
its sum over i is unity, is a measure of how “perpendicular” the i-th species 
axis is to the surface defined by the m-th equation of state in y space. A 
species k is identified as a CSP radical associated the m-th exhausted mode 
whenever Q m (k) is not a small number. 

For the example treated in §3.3, we have for mode # 1 : 


a^ 1 = 


1 


K t +AA 


4 A - 2 K x 0 \ 
-2 A Kx 0 
—2AAHi KxkHx 0 / 


Qi(0 = 


Kx+AA 


The set of its diagonal elements is the radical pointer: 

1 


[4 A, Ki, 0]. 


(65) 


( 66 ) 


This radical pointer informs us that either A or B, but definitely not C, may 
be used as a radical. In other words, only the ODE of either A or B can be 
replaced by (35b). 

Maas and Pope (1992) did not realize that the choice of what species to 
solve for from the equations of state is restricted, and suggested that they 
. . can be chosen quite arbitrarily” with special caution on uniqueness of 
solutions. This is not our experience (Lam and Goussis, 1991b), and (66) 
is developed to provide a quantitative criterion for radicals. We have tested 
our radical pointer criterion by purposely solving for the “wrong” species — 
according to our radical pointer — from (58) in a number of test problems, 
and obtained wrong answers as expected. 



Appendix I 


46 


6.5 The Radical Correction 

What if the values of f™ in (58) are small but not as small as we would like 
at some moment in time? Theoretically, all M modes should rapidly decay 
to some much smaller values in the next several r(Af) seconds, while the 
main reaction activity proceeds with the current characteristic time scale of 
r(M + 1) seconds. If we are not interested in the details of the decay process 
but are only interested in finding an approximate “initial” condition for the 
next time period, we can use Newton’s method to find the change to the 
value of y, Ay rc , which would zero out the values of /™, m = 1,2, . . . , M : 

M 

Ay„ » - £ (67) 

m,n=l 

In other words, y — ► y + Ay rc as the amplitudes of the fast modes decay 
toward zero in the next several r(M) seconds. We shall call Ay rc the radical 
correction. In Lam and Goussis (1991a), the radical correction was referred 
to as “inclusion of the homogeneous solution” in evaluating the change of y, 
but no details were given there. 

When a simplified model is used to compute for an approximation solution 
by neglecting the M exhausted modes, the values of /™’ s are small at the 
beginning of the time period and are theoretically expected to remain small 
in the slow evolutionary period. The radical correction can be used to ensure 
that the numerical solutions adequately satisfy the approximate equations of 
state as required. 

Suppose the initial values in the example worked out in §3.3 does not 
satisfy the partial-equilibrium of mode #1: 

/;(0) = F‘(0) + ^ 0, (68) 

and the solution in the rapid transient period is of no interest to the in- 
vestigator. The radical correction can be used to obtain the effective initial 
values, the counter part of (14), for the next time period. With M = 1, 
t, 1 = + 4v4(0)) -1 , (67) gives: 

2/o 1 (0) 

k x {K x +lA{ 0)Y 


A(0+) « A( 0) - 


(69a) 
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fl(0+) « B{ 0) + 
C(0+) « C(0) + 


PA 0 ) 

k\{Ki + 4A(0))’ 
Ai7i/ O x (0) 
* 1 (A' 1 +4A(0))‘ 


(69b) 

(69c) 


It can easily be verified that (10) is satisfied and A + 2 B is unchanged by 
this radical correction, while the new value of / o 1 (0+) becomes much smaller 
than before. By recursive applications >1(0+) will converge to (14), and 
/*( 0+) will converge to the nearest stable zero. The computation for the 
slow evolutionary time period can now commence with (36) with the new 
initial values. 


6.6 Explosive Modes 

For chemical kinetics, the eigen-values of J are usually real, and mostly 
negative, signifying decaying modes. Occasionally, some modes may have 
positive eigen- values, signifying “ explosive modes ” which are often of interest 
in the study of “ignition” mechanisms (Trevino and Solorio, 1991; Trevino, 
1991). In our example, whenever A < Ki/(2K 2 ) an explosive mode exists 
which is clearly the manifestation of the chain branching of the B radical in 
the reaction system. 

Because the matrix A* is not exactly block diagonal, the non-zero cou- 
pling between the fast exhausted modes and the slow explosive modes can 
cause some exhausted modes to come alive again. From the point of view of 
computation and programming, this complication is easily handled. 


6.7 The Minimum System 

Using the participation and importance indices, it is a simple matter to 
identify the minimum set of species and unknowns (ignore the unimportant 
reactants but includes any species of special interest to the investigator) and 
the minimum set of elementary reactions (ignore the unimportant reactions) 
to form a minimum reduced reaction system which can represent g to any 
reasonable desired user-specified threshold of accuracy for any time period(s) 
of interest. 
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7 Comments on CSP 

In terms of elegance, there is no substitute for a conventional asymptotic 
analysis of an appropriately non-dimensionalized problem with a well de- 
fined small parameter e, with analytical results expressed in terms of simple 
functions. 

In dealing with practical problems, there exists few guidelines for order of 
magnitude estimates and the intelligent non-dimensionalization of the vari- 
ables. The identification of the small parameter e is also not straightforward; 
for a real problem it may be insufficiently small even when successfully identi- 
fied. For moderately complex reaction systems, the formalism of asymptotics 
must be replaced by ad hoc quasi-steady and/or partial-equilibrium hypoth- 
esis, and the role of experience and intuition then becomes central. With 
sufficient algebraic skill, such methodology is indeed capable of generating 
analytical results. However, for massively complex reaction systems, this 
option is simply not viable. 

The theory of CSP provides a new formalism to do asymptotics. Most im- 
portantly, it does not require non-dimensionalization of the variables and the 
identification of a small parameter. The time-dependent value of 
1) is the small parameter. The CSP algorithm can be performed analytical- 
ly for sufficiently simple problems, as was demonstrated in §3.3, and it can 
be programmed to handle massively complex problems — provided the fast 
modes are of the boundary layer type and decay with time eventually. 

CSP can be used to test an intuitive guess, provides a fool-proof way 
to obtain the correct leading approximation in the absence of good ideas, 
and allows the theoretician to concentrate on the task of extracting phys- 
ical insights — because all the massively complicated algebra are left to be 
performed numerically on the computer. 

In most cases, it is relatively easy to interpret a CSP-generated equa- 
tion of state as either a partial-equilibrium or a quasi-steady approximation. 
However, there are also situations when such interpretation is not immediate- 
ly obvious, and the CSP data becomes physically sensible only after further 
reflections. 

From the point of view of asymptotics, CSP is distinctly different from 
the conventional analytical method of “matched asymptotic expansions,” for 
no asymptotic matching is required. A detailed discussion of the differences, 
however, is beyond the scope of this paper. 
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7.1 Local-Eigen Vectors vs. Refined Basis vectors 

As mentioned previously, the set of local eigen-vectors of J is a fool-proof 
choice for use as trial basis vectors. In fact, the special set of basis vectors 
((28a), (28b) and (28c)) used in §3.3 is the leading order approximation 
of the right eigen- vectors of J, with k 2 K 2 / k x Ki playing the role of a small 
parameter. If only so-called leading order approximation (for sufficiently 
small k 2 K 2 lk l Ki) is desired, this set is totally adequate. But what if the 
value of k 2 K 2 /k\K\ of practical interest is only 0.15, and it is desired to 
develop a simplified model with accuracy around 3%? The CSP refinement 
procedure can be recursively applied to improve the accuracy. 

The CSP theory formally guarantees that each full cycle of refinement 
increases the accuracy of the simplified model by the factor t(M)/t(M + 1). 
In the first full cycle, we can generate the leading approximation by using 
constant trial vectors so that no time derivatives are included. In most 
cases, the leading approximation to a complex problem is adequate to provide 
most physical insights. The second full cycle needs the evaluation of first 
time derivatives, and the third full cycle needs the evaluation of second time 
derivatives, etc. Mathematically, these time derivatives are all available. 
Hence CSP is capable of significantly higher accuracy than using just the 
local eigen-vectors. It is, of course, not at all necessary to use the local 
eigen-vectors as trial basis vectors; any informed, intelligent guess can be 
used. In a computer code, it is possible to implement two full refinement 
cycles using relatively little computational resource. 

To continue beyond §3.3, we may choose (28a), (28b) and (28c) as our 
time-dependent trial set, respecting the presence of A(t) in a 2 and b 1 . This 
(second) refinement proceeds routinely except that the time derivative term 
in (41) must be included in evaluating A* . The new CSP-derived simpli- 
fied models so generated will be accurate to 0(k 2 K 2 /kiK\) 2 , compared to 
0(k 2 K 2 /kiKi) for the results obtained earlier. One more cycle of refinement 
will generate models of accuracy 0(k 2 K 2 jk\K x ) 3 , etc. 

7.2 Equation of State and Conservation Law 

In the combustion literature, the concept of conserved scalars is linked with 
the assumption of equal mass diffusivity coefficients and unity Lewis number- 
s. Since no diffusion effects are considered here, these requirements become 
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irrelevant. In this section, we shall show how (approximate) conserved scalars 
for such systems can be found as a by-product of CSP. 

Taking the dot product of b' 0 with (54), we have: 


ffy _ ri 

O © dt /o> 


i = 1,2 


(70) 


If a. row vector b" can be expressed as the gradient (with respect to y) of a 
scalar function 0 n ,o(y) times a scalar function 0 n ,o(y)> then the mode shall 
be called an holonomic mode — a concept well known in classical mechanics. 
The corresponding equation is (70) is called an holonomic constraint, and 
can be rewritten as: 

&n,o ~~TT~ = selected n’s. (71) 

at 


We shall call 0„, o , when it exists, an eligible conserved scalar. 

When the m-th fast mode becomes exhausted, an equation of state is 
obtained by setting f™ & 0. If this mode is also holonomic, the associated 
eligible conserved scalar 0 m , o then becomes a exhausted conserved scalar — 
its value remains constant from the moment of exhaustion. It is likely that 
all exhausted refined fast modes are holonomic, and that 0 m>o and are 
linearly dependent — thus no new information is provided by the exhausted 
conserved scalars. 

When the amplitude of the if-th slow mode is found “negligible,” setting 
this dormant mode fjf « 0 yields no information at all. We can only conclude 
from (70) that the component of the vector dy/dt “in the the direction of 
a ° K r ‘ is negligible. However, if this dormant mode is holonomic (e. g. are 
constants, yielding Bk,o — 1? 0k> = b^ Oy), then the associated Ok,o is a 

dormant conserved scalar— its value is approximately conserved whenever the 
K - th mode is dormant. Conserved scalars in the combustion literature are 
dormant ones, and they almost always have simple physical meanings, such 
as conservation of atomic species or total energy. While &k,o ~ constant is 
also an algebraic relation between the species, we do not call it an equation 
of state because the constant is not universal. 

In our example, the basis vectors used in §3.3 are the leading order results 
obtained by refining the trial basis vectors given in §3.2. Equations (33a), 
(33b) and (33c) correspond to (70) and they are all holonomic. The eligible 
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conserved scalars are (M = 1): 

01.0 = -A 2 + KiB, (72a) 

02.0 = A + 2B, (72b) 

© 3i o = (AH 1 -AH 2 )A + (AH l -2AH 2 )B + C, (72c) 


^l.o — ^2,o — 1> ^3,o 1 4 

K i -f 4A 


All three refined modes in this example are exactly holonomic. In general, 
there is no guarantee that any CSP refined mode is holonomic 4 — it is not 
possible to determine from numerical CSP data whether a mode is holonomic 
or to deduce the form of the possible eligible conserved scalar except for the 
obvious case of time-independent b". 


7.3 The User-Specified y error 

The role of a theory is to construct a model which should be as simple as pos- 
sible and yet can generate valid predictions in some domain of the parameter 
space of interest. A good model need not predict all things correctly. For 
example, a certain simplified model may be very good indeed for the heat 
release and temperature history, but very bad for the predictions of certain 
pollutant concentrations. For the same reaction system, different simplified 
models may be appropriate depending on the desired scope and accuracy of 
the solutions. 

The CSP theory addresses this issue by choosing not to directly define a 
vector norm for the A- dimensional state space of y. Instead, it asks the in- 
vestigator to supply the vector y error , which specifies the amount of absolute 
error considered tolerable for each component of y. This step is intuitively 
sensible, and is conceptually equivalent to choosing a norm. The CSP the- 
ory relies on this y error to decide whether a decaying fast mode should be 
neglected. As a consequence, a CSP-derived simplified model which focuses 

4 A sufficient but not necessary condition for the n-th mode to be holonomic is that 
b n = constant. Hence contant trial row basis vectors always generate holonomic modes, 
but the resulting eligible conserved scalars may never become conserved scalars. For our 
example, if the default trial b 1 in §3.2 were used, we would obtain 0i = —A — B which, 
in contrast to 0j (O , is not conserved after fast mode #1 is exhausted. 
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on temperature only may be quite different from one which claims additional 
predictive capability such as pollutant concentrations. 

If too stringent a y err or is chosen, CSP would simply report M = 0 
and would fail to recommend any simplified model. As y erTOT becomes more 
forgiving, CSP would recommend simplier models consistent with the spec- 
ified tolerance. This posture of producing simplified models on demand is 
refreshing — when compared to the posture of conventional asymptotics. 


7.4 Insights from CSP Numbers 

In a theoretical analysis, no specific numerical values of the parameters in 
the problem need be specified; only their order of magnitude is assumed. 
While the accuracy of the analysis may be limited — because the small pa- 
rameter may not be that small in practice — the parametric dependence of 
the results are available for inspection and interpretation. An experienced 
theoretician can speak knowledgeably about the reaction system and inter- 
pret the analytical results in physically meaningful and intuitive terms such 
as chain branching, radical formation, ignition temperature, alternative re- 
action paths, etc. 

When results are obtained computationally, each set of numerical data is 
valid only for one specific set of parameters. For massively complex problems, 
the computer printouts are overwhelming. The interpretation of these num- 
bers is a much more difficult matter. The CSP data generated, which is all 
numerical, is designed to be easy to use to extract qualitative and physically 
meaningful insights. 

For non-linear problems, the CSP data is itself massive. However, it need 
not be computed at every time step; it is needed only when the theoretician 
wants to know: what is going on here? While the time dependence of the 
reaction system may appear to be very complex, the CSP-derived refined 
basis vectors may depend only weakly on time, and the simplified model 
derived at selected points in time can provide insights about the “physical 
mechanism” of the system for a finite time interval. In fact, significant change 
of the behaviors of basis vectors is a signal that the physical mechanism has 
changed, and a new simplified model is needed. For linear problems, a single 
set of CSP data — the eigen structure of the constant matrix J — will explain 
everything for all time. For our simple example (which is non-linear), a single 
set of time-independent trial basis vector is adequate to generate the correct 
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refined set for all time. 

From the CSP viewpoint, the problem of simplified kinetics modeling re- 
duces to finding a set of basis vectors which make A‘ block-diagonal. The 
conventional methods guess at them, while CSP provides a rational and pro- 
grammable way of refining a trial set to make the resulting A) more diagonal 
than before. For massively complex problems, the computer is assigned the 
task of computing an adequate set, and the investigator uses common sense 
to extract physically interesting information from the computer generated 
numerical basis vectors. The participation index and the importance index 
are provided to quantify the roles played by each elementary reaction includ- 
ed in the model. Together they provide answers to most of the interesting 
questions about a reaction system. Any simplified model derived based on 
intuition and experience can now be compared with the CSP-derived results 
to verify that the various judgments and guesses used are correct. 

Depending on the user-specified accuracy threshold, y error , certain chem- 
ical species and elementary reactions can be removed from the reaction sys- 
tem, producing the so-called reduced reaction system. The more lenient y err or 
is, the simpler the CSP-derived reduced system — while the essential charac- 
ter of the full system is retained. In Goussis and Lam (1992), the methanol 
oxidation problem was studied with 30 species plus temperature and 173 el- 
ementary reactions. The CSP data generated can be roughly divided into 
four epochs, each characterized by a distinctly different set of CSP-derived 
basis vectors, and together they indicated (through the CSP indices) that 
a reduced reaction system with 15 species (including methanol) would have 
adequate accuracy for all 15 unknowns. This CSP-derived insight was com- 
putationally confirmed. This reduced reaction system, of course, makes no 
predictions on the discarded species. 

7.5 CSP vs. Sensitivity Analysis 

The standard method of sensitivity analysis (Yetter, Dryer and Rabitz, 1991) 
is to linearize the problem and then to compute for the linearized response 
to various perturbations. For example, if one wishes to know what happens 
to y 95 when k 12 7 is changed by a small amount, one numerically computes 
d(\ogy 9S )/d{\ogk n7 ). 

CSP provides an alternative to this brute-force method. First, we look 
at the participation indices Pj 5 7o , m = 1,2, , M. If these indices indicate 
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that reaction #127 is a significant participant of the m-th exhausted fast 
mode, then k 127 will directly affect the values of the radicals pointed to 
by the radical pointer Q m (i)- Otherwise km has no effect on the radicals. 
Next, we look at the importance index I^m', its numerical value tells us how 
strongly km affects the reaction rate of j/ 95 . 

Frequently, it is desired to know what would happen if additional ele- 
mentary reactions were added to a reaction system after a computation is 
already done. CSP can be used to provide useful answers. To illustrate, let 
reaction #4 be added to our example reaction system: 

reaction # 4 : B + B ^ A. (74) 

Its stoichiometric vector and reaction rate are: 

s 4 = [l, -2, A tf 4 ] r , F 4 = k 4 (B 2 — K 4 A). (75) 

Since s 4 is not linearly independent, we have: 

s 4 = Si - 3s 2 - (A Hi - 3A H 2 - AH 4 )s 3 . (76) 

Hence: 


B\ = B[ - ZB\ - (A#i - 3A#2 - A/f 4 )J3^, (77) 

<t w = - Kf w - (Atfj - 3Atf 2 - AH 4 )s° 0 f' \ (78) 

Using the basis vectors in §3.3, we can computed the following CSP data: 



d 1 _ 1 di — Ki±M 
D \ ~ 1 ’ °2 Ki+iA’’ 

£] = o, 

(79a) 


B\ = 0, B\ = 1, 

Bl = 0, 

(79b) 


B\ = 0, B\ = 0, 

Bl = 1, 

(79c) 

and (for 

the time period when M = 1 only), 



OySloW 
S o, 1 

II 

O 

o 

,o 


(80a) 


s °2 OW = 2A ’ (Ki + 4A)AH 2 -(K 1 +2A)AH 1 ] t , (80b) 

Ki + 4/i 

s° 0 % ,ow = [0, 0, 1] T . (80c) 
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The effective stoichiometric vectors for the other two time periods are also 
available. It is now a simple matter to compute the participation and im- 
portance indices for reaction #4 from this CSP data and decide whether its 
previous omission can be justified in any time period of interest. If the in- 
dices P\ 0 and f 4 ° show that reaction #4 is a major player, then the old basis 
vectors should be refined to reflect the new character of the new reaction 
system. 

7.6 Effects on Diffusion 

If spatial diffusion is included, the governing equations will contain the Lapla- 
cian operator, and become a system of PDE’s: 

| = g + D0V 2 y (81) 

where D is — in most cases — a N x N diagonal matrix of diffusion coefficients 
which are considered constants here for the sake of simplicity. 

Assuming that the spatial diffusion mechanism is not the fastest process, 
we can obtain the corresponding approximate equations of state, (58), when 
the fastest M modes are exhausted as before. In the slow evolutionary time 
period and spatial domain, the fast modes can be neglected from the right- 
hand side of (81) to obtain, instead of (60), the following simplified model: 

^ w £ s ° 0 f w F r + D: O V 2 y (82) 

r=l 

where D° is the effective diffusion matrix , the projection of D in the slow 
subspace, given by: 

M 

d; = (I- £<Cbr)0D- (83) 

m=l 

Let the diffusion matrix D in our example be diagonal, with Da, Db, 
and Dc as the diagonal elements: 

l D a 0 0 \ 

D = 0 D b 0 . (84) 

V 0 0 D c ) 
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When M = 1, the effective diffusion matrix D° can be computed from (83) 
and (65): 


D° = 

0 K x + 4A 


K x D a 

2AD a 

2AAH x D a 


2K x D b 

4AD b 

-KrAHxDg 


0 

0 

{K x + 4A)D C 


(85) 


This derivation of D° using CSP basis vectors is straightforward and pro- 
grammable, and can be used to deal with massively complex PDE systems. 
Adjacent to solid boundaries, spatial boundary layers will exist, and effective 
boundary values for (82) are needed and must be carefully derived. 


8 Discussion 

Experience and intuition have always been central to the conventional deriva- 
tion of simplified models. It was generally known that both the partial- 
equilibrium and the quasi-steady approximations were related to convention- 
al analytical singular perturbation procedures — provided that the relevant 
small parameter(s) could be identified. 

The goal of developing a general theory of singular perturbation which 
can handle any system of non-linear first order ODE’s in a programmable 
manner without the identification of small parameter(s) seems very ambitious 
indeed. The CSP theory developed by Lam et. al. has achieved this goal, 
but only for boundary-layer type problems where all fast modes eventually 
decay exponentially. Fortunately, most problems in chemical kinetics are of 
this type. 

The basic strategy of CSP is to uncouple the fast, exhausted modes from 
the slower, currently active modes though an intelligent choice of basis vec- 
tors. For a strictly linear problem, the ideal uncoupling basis vectors are the 
time-independent eigen-vectors of J. For non-linear problems, the desired 
uncoupling basis vectors are time-dependent because A } depends on y (t). 
The main concession made by CSP is to abandon the goal of diagonalizing 
Aj and be satisfied with a nearly block- diagonal one. A programmable re- 
cursive refinement procedure is provided to successively weaken the fast-slow 
couplings. 

Pragmatically, the extraction of physical insights from the CSP data is 
advantageous only if the CSP data, such as the participation index and the 
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importance index, changes less rapidly with time than the solution y(t). In 
problems involving chemical reactions, the basic chemical process is usual- 
ly described by a reaction rate F* and a constant stoichiometric vector s r . 
Whatever simplifications which are available can usually be attributed to 
some reactions being much faster than others and had now spent themselves. 
As a consequence, good decoupling basis vectors in chemical kinetics are 
most likely to be weakly time-dependent, because the relevant chemically 
meaningful stoichiometric vectors of the fast reactions are time-independent. 
In the example studied in §3.3, the information extracted from the CSP data 
hardly change with time at all, because the “mechanism” of the problem did 
not change. 

The CSP algorithm depend on the fast modes to become exhausted, which 
translates into assuming that the fast eigen- values of J are “essentially real” 
and negative. Theoretically, if a highly oscillatory fast mode does not damp 
out (i.e. a Wentzel, Kramers and Brillouin (WKB) type problem), CSP 
will only be able to identify this mode, but will not be able to provide an 
algorithm to handle the evolution of y in the next larger time scale. The 
extension of CSP to WKB type problems remains to be explored. 

CSP can be used, in principle, on problems which have a large number 
of modes but only moderate amount of time scale separation (i.e. the ratio, 
(t(M + 1 )/r(M)), is only moderately large), provided that only moderate 
accuracy is desired. It thrives on large time scale separations, and performs 
best when the separations are “asymptotically large.” 


9 Concluding Remarks 

A computer code, called CSP91 and programmed by Dr. Dimitris Goussis, 
has implemented the CSP algorithm, and is available on request. The code 
can be used either as a diagnostics subroutine, or as an stiff ODE solver. It 
has been used to demonstrate in a number of test problems (Lam et. al., 1989; 
Goussis et. al. 1990; Goussis and Lam, 1992) that numerical solutions gen- 
erated by the CSP-derived simplified models are in quantitative agreement 
with that of the corresponding full models in accordance to the accuracy 
thresholds specified. Different simplified models are generated for different 
user-specified y „ T0T . The CSP-derived insights, such as which fast reactions 
are in partial-equilibrium with each other (through the participation indices), 
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which reactions are rate-controlling (through the importance indices), which 
species can be considered as radicals in the CSP context (through the rad- 
ical pointers) and how they can be solved for from the equations of state 
(through the radical correction), have been thoroughly tested and verified. 
In most cases, the CSP-derived insights are consistent with the expectations 
of competent chemical kineticists knowledgeable in the reaction system under 
study. For the few cases when the CSP-derived insights appeared surprising 
at first sight, they all eventually became obvious after further reflections. 
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A STUDY OF HOMOGENEOUS METHANOL OXIDATION KINETICS USING CSP 

D. A. Goussis and S. H. Lam 
Princeton University, Princeton, NJ 08544 


Abstract 

The homogeneous oxidation of methanol in air at constant pressure is 
examined using data generated by the method of computational singular perturbation 
(CSP). At any moment in time, the number of exhausted fast modes and the radicals 
(sometimes called the intermediaries) are computationally identified. The 
participation index, which quantifies the participation of any elementary reaction to 
an equations of state of the radicals, along with the importance index, which 
quantifies the importance of any elementary reaction to a particular species of 
interest, are computed and used to assess the sensitivities of the solution to the 
reaction rate constants. Every elementary reaction is classified so that it either 
belongs to the equilibrated set which contains fast reactions already equilibrated 
among themselves, and/or the rate-controlling set which contains reactions 
controlling the current rate of activities, or neither of the above sets - in which case 
it is superfluous. A number of numerical experiments were performed to verify the 
assessments: (a) the relative effectiveness of the reaction rate constants of two 
reactions (#16, #160) in breaking up the fuel indicated by the importance index is 
verified, (b) that fuel breakup in an early time period can actually be slowed down 
by increasing the reaction rate constants of certain fuel breakup reactions (#156, 
#159) is verified. Numerical experiments also show that species identified as 
radicals respond instantly to sudden changes in reaction rates, while the non- 
radicals respond more smoothly. The overall response of the unknowns to 
perturbations is always consistent with the CSP-derived effective stoichiometric 
coefficients. In addition, a minimum set of species is constructed with the help of 
the CSP data. This minimum set, which trims the original full set of 30 species to 
15 species, generates numerical solutions in excellent agreement with solutions 
obtained with the full set 


1. Introduction. 


The present paper studies the constant pressure oxidation of methanol in air. The full 
kinetics mechanism, taken from Egolfopoulos, Du and Law 1 [EDL], consists of 30 species and 
173 reversible elementary reactions. A partial list of the elementary reactions is given in Appendix 
I. The case of a fuel-lean (equivalence ratio=0.6) mixture at 1 atmosphere will be considered. The 
initial state of the system is taken to be: T(0) = 1027°K, Y[CH 3 OH](0) = 0.00779, Y[O 2 ](0) = 
0.01980, Y[N 2 ](0) = 0.9724, where T is temperature and Y is mass fraction. The calculation is 
performed in a Chemkin environment 2 , and the resulting numerical solution is analyzed using the 
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data generated by the method of computational singular perturbation (CSP) developed by the 
authors 3 ’ 4 * 5 . 


Figs, la and lb show the time evolution of four species; similar plots for all the other 
species and the temperature are readily available. The oxidation of CH 3 OH can be described as a 
four-stage process. The first stage, 0.00 seconds < t < 0.05 seconds, is an incubation period in 
which certain radicals are created. The breakup of the fuel occurs in the second stage, 0.050 
seconds < t < 0.079 seconds. The third stage, 0.079 < t < 0.081, consists of the very rapid 
conversion of CO to CO 2 . The conversion of the remaining CO to CO 2 takes place in the fourth 
stage, 0.081 < t < 0.130, with progressively slower rates. The integration of the kinetics equations 
is straightforwardly performed by a CSP ODE solver which, in addition to generating solutions of 
guaranteed user-specified accuracy, also generates a set of CSP data. Each time interval between 
markers shown in Figs. la,b covers 24 integration time steps selected by the CSP code. 

2. Preliminary Discussions 

The conventional derivation of a simplified model requires that the "radicals" and fast 
reactions in the reaction system be somehow identified, and the quasi-steady and partial- 
equilibrium approximations be systematically 6 ’ 7 ’ 8 ’ 9 applied. Considerable intuition, experience 
and mathematical skills are required. The simplified model so obtained is highly valued because it 
identifies the rate-controlling elementary reactions, and can be used to provide or interpret the 
sensitivity information 10 ’ 11 ’ 12 of the reaction system. The CSP method performs the above tasks 
routinely using a programmable computational algorithm, and can generates time-resolved 
simplified models without the need of intuition and experience. A brief summary of the CSP 
method is given in Appendix D. 

3. The Identification of Radicals 

Including temperature, the system of ODEs of the [EDL] full mechanism consists of 3 1 
unknowns and 173 reversible elementary reactions. Let y = [Y 1 , Y 2 , ... Y N ] T be a column vector 
of the N = 31 unknowns. Let F^(y) £ 0 and F k (y) £ 0 be the forward and backward reaction rates 
of the k-th elementary reactions, where k=l, 2, ... R with R=173. In general, the governing 
equations can be written is the following compact form: 


( 2 . 1 ) 
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where s, is called the stoichiometric vector of the k-th elementary reaction, and the summation 
£ 

convention is used. Since this is a 31 -dimensional system, the solution is expected to have 31 
"modes.” These modes are ordered by CSP according to their speed: mode #1 being the fastest 
and mode #26 being the slowest. Modes #27 to #31 are constants and do not evolve with time; 
they represent the conservation of energy, the conservation of C, H, and O atoms and the inertness 
of N 2 . In addition, CSP computationally determines M, the number of exhausted fast modes (see 
Appendix II). Fig. 2 shows M vs. time. Each of the exhausted fast modes is associated by CSP 
with either a quasi-steady approximation (for a certain radical species) or a partial-equilibrium 
approximation (for a certain set of fast reactions). In the present paper, we shall call a species a 
radical when it can accurately be computed in terms of the other "major species" from certain 
approximate "equations of state" derived by CSP. More will be said on this in §4 later. The 
identification of these radicals was achieved by the use of Q(m; i), the "radical pointer" of each of 
the m exhausted modes, as described in Appendix II. 

The following information is obtained from the CSP data for each of the four stages: 

First Stage to 00 < t < fi 05 seconds'). M = 14 . The following 14 species are identified as radicals: 

ch 3 , c 2 h 5 , ch 3 co, hco, h, ch 2 , c 2 h 3 , ch 2 oh, hcco, o, ch, c 2 h, ch 3 o, oh. 

Second Stage (fi 050 < t < 0.079 seconds'). M = 15 . HO 2 joins the above, making a total of 15 
radicals. The breakup of the fuel is in full swing. 

Third Stage to 079 < t < 0,081 seconds'). M = 12 . This is a most active period in which things 
happen very rapidly. Three exhausted modes becomes alive again, and CH 3 , H and O are 
temporarily removed from the list of radicals. 

Fourth Stage. (0.081 < t < 0.130 seconds'). M = 15-24. After the frenzy activities in the third stage, 
CH 3 , H and O rejoins the list of radicals. As the fourth stage progresses, M increases from 15 to 
24, mostly as the results of extinction of the carbon related species (CO and CO 2 excluded). The 
remaining non-carbon related radicals in this late period are O, H, OH, HO 2 and H 2 O 2 . 

4. The Simplified Model af t=0.0345 

We shall demonstrate in some detail the procedures for extracting information about the 
reaction system using the CSP data in one typical moment in time. At t = 0.0345, when the fuel is 
breaking up and HO 2 is not yet a radical but is about to peak, CSP determined M = 14, and 
generated 14 approximate "equations of state": 
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m =1, ... , 14. (4.1) 

Each of these equations of state can be physically interpreted either as a quasi-steady approxi- 
mation or a partial-equilibrium approximation. Samples of these 14 equations of state, suitably 
normalized, are displayed below: 


Exhausted Mode #1 [CH 3 ]: 
F 29 


f 33 + f 32 + f | 63 + 


0.50 » 0.43 +0.05 +0.01 + ... , 


(4.2a) 


Exhausted Mode #5 [H]: 


-43 -53 _3 
f + +f; + F + 


fJ . 56 + f \ 2 + F + 8 + f | 57 + fJ. + f ! 1 +F 9 + 


0.35+0.14 +0.02 - 0.23 40.08 +0.08 +0.06 


+ ’ * + 
40.02+0.02 +0.02 + 


(4.2b) 


Exhausted M ode #13 [CH 3 O]: 

F|6° + Fi 57 . F 43 + 

0.42 +0.08 * 0.48 + ... . 


(4.2c) 


Each exhausted mode is identified by the radical pointer Q(m; i) described in Appendix II with a 
specific species which shall be referred to as a radical. The left and right hand sides of (4.2a, b,c) 
are the rates of production and consumption of the identified radical. The terms on each side are 
ordered in descending order of magnitude; sufficient terms are kept so that the total contribution of 
the omitted terms is below a user-specified threshold. The numbers displayed below each term is 
the participation index, P™, which is defined in Appendix II and measures the significance of the 
k-th reaction to the m-th exhausted fast mode. 


In addition, CSP derives the following simplified model: 

f-cX-F*). <«> 

where is called the effective stoichiometric vector of the k-th reaction. Mathematically, it is the 
projection of the original stoichiometric vector in the currently active slow subspace. Samples 
of these 31 equations are displayed below: 


d[CH4] 
dt 


17^16^ 17 _160 
+ c 160 F + 


" c 16 F + 


17^9 


+ C29P. + c^F^ + cJ5 6 F 


.43 

'+ 


17 


156 

+ 


17-33 , 17 -159 
+c 33 F + +c 159 F + 
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+ 




[CH4]=Y17, 


(4.4a) 


d[H] 
dt 


7 -16 . 7 -15 


= c 16 F + 


+ c 15 F + 


+ c 162 F + 


162 


[H]=Y?, 


(4.4b) 


d[CH 3 OH] 
dt 


15-16 15 -160 15-40 

+ c 160 F + + c 40 F + 


= c 16 F + 


0.40 - 0.16 + 0.12 


15 -162 
+ c 162 F + +c 
-0.05 


15-57^ 15-53 
57* + +c 53 F + 

+ 0.05 + 0.05 


+ ... , [CH 3 OH] = Y15, (4.4c) 

where the coefficients c\. (the i-th element of c fc ) are provided by CSP numerically. The numbers 
displayed below each term in (4.4c) is the importance index, 1^ , a measure of the significance of 
the k-th reaction to the simplified equation for CH3OH. The definition of 1^ is given in Appendix 
n. Again, the terms on the right hand side are ordered in descending magnitude; sufficient terms 
are kept so that the total of the omitted terms is below a user-specified threshold. Note that the 
solution of these 31 equations automatically satisfies the 14 approximate equations of state given 
symbolically by (4.1) and explicitly by (4.2). Once may freely replace the differential equation for 
any radical, such as (4.4b), by its equation of state, in this case (4.2b). 


It is interesting to compare the CSP derived simplified model equations with the original 
equations. For [CH 3 OH], the original equation is: 

32,-159 
P 

where p is total mass density. Using the CSP-computed numerical values for the simplified 
equation for [CH 3 OH] is: 


( ! 


drCH 3 OHh 

dt Joriginal 


,160 
■4 

+ 


.156 
+ 


- — {f;~ + + Fr~ + Fr~ + Fr + 


.162 
4 

+ 


.157 
*+ 




(4.5a) 


,15 


^d[CH 3 OHl^sp . . — {2.25F^ + 1.19F+ 60 - 0.75F* 0 + f | 62 - 0.34FJ 7 + 0.94F^ 3 


+ ...). (4.5b) 

The terms in (4.5a, b) are ordered in descending magnitude. In contrast to (4.5a) where each term 
involves methanol as a participant in a reaction, only reactions #160 and #162 in (4.5b) directly 
involve methanol. Similar observation can be made for all the simplified equations. For example, 
none of the seven reactions listed in the CSP-derived (4.4a) for [CH4] involves CH4. Theoret- 
ically, the difference between the right hand side of the original and the CSP-derived simplified 
equations is some linear combinations of the equations of state given symbolically by (4.1) and 
explicitly by (4.2). Even though the CSP-derived simplified equations frequently do not make 
chemical sense at first glance, our experience is that they usually make sense upon further 
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reflections. Computationally, solutions obtained from the CSP-derived simplified equations and 
from the full equations agree within user-specified bounds. 

Similar information can be obtained for any moment in time. 

5. The Importance Index 

The importance index 1^ defined in Appendix II measures the importance of the k-th 
reaction to the i-th unknown. Figs. 3a and 3b show I^ 5 , the importance index for CH 3 OH, for 
reactions which do and do not involve CH 3 OH, respectively. These two graphs provide concise 
information on the rate-controlling reactions for the breakup of CH 3 OH. By inspection, Fig. 3a,b 
show that reaction #16 is the most effective reaction in both the first and the second stages, even 
though it does not directly involve CH 3 OH. It is seen to be somewhat more effective than 
reactions #160 and #162 which do. A numerical experiment was performed by increasing the rate 
constants for reactions #16 and #160 by a factor of 2.5 in the first stage. The computed time 
histories of CH 3 OH confirmed the above qualitative prediction. 

In addition. Fig. 3a shows that reactions #156 and #159, which chemically consume 
CH 3 OH, have positive importance indices in the 0.00 < t < 0.03 period: ljj 6 > 0 and ij^ > 0. In 
other words, they effectively contribute positive terms in the first stage to the right hand side of 
(4.4c) - a somewhat counter-intuitive prediction. Another numerical experiment was performed in 
which the rate constants for reactions #156 and #159 were simultaneously doubled. The resulting 
data shows that the breakup of CH 3 OH was indeed slowed and delayed in the first stage, again 
confirming the CSP expectation. In the third stage, however, faster #156 and #159 do favor the 
breakup of CH 3 OH, as indicated also by Fig. 3a. 

6. Classification of Reactions 

The left hand side of (4.1) consists of positive and negative terms which nearly cancel each 
other. Separating the positive and negative terms and placing them across an equal sign as is done 
in ( 4 . 2 ) for each of the exhausted modes, we can pick out reactions which participate most 
significantly on each side using the participation index, P™, defined and described in Appendix II. 
We shall call the total collection of such reactions the equilibrated set . Loosely speaking, these are 
the fast reactions which have equilibrated among themselves. At t = 0.0345, the equilibrated set 
consists of: forward 1, 3, 9, 11-14, 16, 30, 32, 33, 35, 38-41, 43, 53, 57, 65, 67, 71-74, 78, 
84, 86 , 94 , 98, 101, 109, 110, 116, 120, 123, 128, 132, 141-143, and backward 29, 113, 118, 
142, 147, 149. 
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We can pick out, for each unknown, reactions which contribute most significantly to the 
right hand side of (4.3) using the importance index, 1^. We shall call these the rate-controlling set. 
Loosely speaking, these are the slower reactions currently driving the system. At t = 0.0345, the 
rate-controlling set consists of: forward 12, 14-16, 33, 38, 40, 41, 43, 53, 57, 58, 60, 62, 101, 
109, 110, 127, 156, 157, 159, 160, 162, and backward 29, 102, 127, 136. Note that these two 
sets are not totally distinct; some reactions ( e.g . 12, 14, 16, etc.) belong to both sets. A reaction 
which belongs to neither set is obviously not important to the system and is therefore superfluous 
at that moment in time. 

The reactions which only belong to the equilibrated set define via (4.1) the so-called slow 
subspace which may be visualized as a (N-M)-dimensional "surface" on which the N-dimensional 
solution point y(t) moves. Changing the rate constants in the equilibrated set changes the surface, 
but does not change the speed of the motion. The reactions which only belong to the rate- 
controlling set control the speed of the motion of y(t), but do not affect the surface on which it 
moves. The reactions which belong to both sets affect both. Numerical experiments were 
performed by perturbing the reaction rate of reactions #15 in a step-function manner for t ^ 
0.0345. Since #15 belongs only to the rate-controlling set, the response of y(t) to this abrupt 
perturbation is smooth and agree qualitatively with the CSP expectations, additional numerical 
experiments were performed by perturbing the reaction rate of reaction #160 which belongs to 
both sets. The response of y(t) is again in qualitative agreement with the CSP expectations: some 
radicals respond discontinuously, all major species respond smoothly. Such information can be 
very useful in reaction-path analysis. 

7. The Minimum Set of Species 

Normally, one is usually interested only in the time history of a few species in a reaction 
system. Using the data generated by CSP, it is straight-forward to use the CSP data to identify the 
minimum set of unknowns (species) which includes an user-specified species of interest. The 
following minimum set which includes CH 3 OH is obtained for the time interval studied: CH 3 OH, 
CH 2 0, H, H 2 0 2 , CH 3 0, CH 3 , H 2 , H0 2 , CH 2 OH, CO, O, H 2 0, HCO, N 2 , 0 2 , OH. The 
solution of this minimum set of 16 species (plus temperature) system is in excellent agreement 
with the solution computed from the full reaction system. Comparison of results obtained for CO 
and H are shown in Fig. lb. 
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8. Conclusions 

Neither of the authors are knowledgeable about chemical kinetics in general, and about 
methanol oxidation in particular. The above information are presumed by us to be informative, and 
are obtained routinely using the CSP data generated by the CSP code. Because of limitation of 
space, only part of the available information is included here^. The method is clearly useful for 
the study of massively complex systems when intuition and experience are lacking and 
conventional analysis is untenable. 
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Appendix I 

Selected Elementary Reactions From Full Kinetics Mechanism of [EDL] 


(I) h + o 2 = o + oh 

(3) OH + F^ = H + I^O 

(4) OH + OH = O + F^O 
(9) H + 0 2 + M = H0 2 + M 

(II) ho 2 + h = o 2 + H2 

(12) H0 2 + H = 0H + 0H 

(13) H0 2 + 0 = OH + 0 2 

(14) H0 2 + OH = 1^0 + 0 2 

(15) H0 2 + H0 2 = W 2 0 2 + 0 2 

(16) H 2 0 2 + M = 20H + M 

(29) CH 3 + OH = CH 3 0 + H 

(30) CH 3 + OH = CH 2 + HjO 

(32) CH 3 + 0 2 = CH 2 0 + 0 

(33) CH 3 + H0 2 = CH 3 0 + OH 

(35) 2CH 3 = C 2 H 5 + H 

(36) 20^ (+M) = C 2 H 6 (+M) 

(38) CH 2 0 + H = HCO + ^ 

(39) 0^0 + O = HCO + OH 

(40) ch 2 o + oh=hco + h 2 o 

(41) CH 2 0 + H0 2 = HCO + Kp 2 
(43) CH 3 0 + M = CH 2 0 + H + M 
(53) HCO + M = CO + H + M 

(57) HCO + 0 2 = CO + H0 2 

(58) HC0 + H0 2 = C0 2 + 0H + H 
(60) C0 + 0H = C0 2 + H 

(62) CO + H0 2 = C0 2 + OH 
(65) CH + 0 2 = HC0 + 0 
(67) CH 2 + H = CH + H 2 

(71) CY^ + OH = CH + H 2 0 

(72) CH 2 + 0 2 = C0 2 + H + H 

(73) CH 2 + 0 2 = CT^O + O 

(74) CH 2 + 0 2 = C0 2 + ^ 

(78) C^ + Kp 2 = CH 3 O + OH 


(84) C 2 H 6 + H = C 2 H 5 + H 2 
(86) C 2 H 6 + OH = C 2 H 5 + H 2 0 
(94) C 2 H 5 + 0 2 = C 2 H 4 + H0 2 
(98) C 2 H 4 + H = C 2 H 3 + H 2 

(101) c 2 h 4 + OH = c 2 h 3 + h 2 o 

(102) C 2 H 4 + OH = CHg + 0^0 

(109) C 2 H 3 + 0 2 = C 2 H 2 + H0 2 

(1 10) C^ + 0 2 = 0^0 + HCO 
(113) C 2 H 2 + H = C 2 H + ^ 

(i 16) c 2 h 2 + oh = c 2 h + up 

(120) C 2 H + 0 2 = CO + HCO 
(123) CI^CO + H = HCCO + ^ 

(127) CH 2 CO + OH = CH 2 0 + HCO 

(128) CH 2 CO + OH = HCCO + F^O 
(132) HCCO + 0 2 = 2CO + OH 
(136) CHjHCO = CHj + HCO 

(142) CI^CO + M = Cl^ + CO + M 

(143) Cl^CO + H = CHg + HCO 

(147) Cl^CO + = Cl^HCO + H 

(154) O^OH (+M) = CI^ + OH (+M) 

(155) CF^OH (+M) = CHjOH + H (+M) 

(156) CH 3 OH + H = CHjOH + F^ 

(157) CH 3 OH + H = CH 3 0 + ^ 

(158) CH 3 OH + O = CF^OH + OH 

(159) Cl^OH + OH = CF^OH + F^O 

(160) CH 3 OH + OH = 0^0 + I^O 

(161) CHjOH + 0 2 = CF^OH + H0 2 

(162) CHjOH + H0 2 = CH 2 OH + YLp 2 

(163) CH 3 OH + CH 3 = CF^OH + CH 4 

(164) CI^OH + CH 3 = CHjO + CH 4 
(167) CHjOH + 0 2 = 0^0 + H0 2 
(170) CH 2 OH + CHj = C 2 H 5 + OH 
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Appendix II 

A Summary of CSP Concepts 
A1 Basic CSP Concepts and Data 

The basic CSP idea is to split the N-dimensional space of the vector g into two subspaces, 
a fast and a slow subspace: 

g - gfast + gslow^ ( A 1 . 1 ) 

where g fast is spanned by a set of M linearly independent column basis vectors a m » (m=l, 2, ... 
M). CSP provides an algorithm 5 to determine M and to compute for this set of a along the 

jH # # 

trajectory of the solution y(t), together with a set of row vectors b which is ortho-normal to a m : 

b m • a = 6™, m, n = 1,2, M. 
n n 

The time-resolved values of M, the basis vectors a m and b m appropriately ordered in 
ascending time scales, are the basic CSP data. 

A.2 Using the Basis Vectors 

The fast subspace, being M-dimensional, contains M fast reaction modes, or simply 
modes. For chemical kinetics problems, these modes are usually decaying modes; i.e. they tend to 
become exhausted. When g fast falls below some user-specified threshold for t ^ tM. we have: 

gfast » o, (A2.1) 

which yields M algebraic relations between the elements of y. In other words, M equations of 
state, usually derived analytically using conventional singular perturbation technique, are directly 
obtained from (A2.1). Another consequence is that an approximate time evolution equation for y 
is obtained for t £ tM- 

& = gslow, (A2.2) 

which is the desired simplified model of the reaction system. CSP guarantees 5 that solutions of 
(A2.2) satisfying (A2.1) at t = tM will satisfy (A2.1) for t £ tM, under fairly general conditions. 

With M and both a and b m available, we can express g fast in terms of them: 
m 

gfast = a b m • g = a f”\ (A2.3) 

m m 

where f 1 " = b m • g = B™ {F^. - F^ } , B™ = b m • s fc . Using (A2.3) in (Al.l) and solving for 
gsl° w , we have: 
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g stow = c k <F$;-F k ), 
c k =a-» m b m )-s k . 


(A2.4a) 

(A2.4b) 


where is called the effective stoichiometric vector, and elements of this column vector, c k , are 
called the effective stoichiometric coefficients of the k-th reaction. 

A.3 Participation Index, Importance Index and the Radical Pointer 


CSP requires the user to specify Ay er ror, the absolute or relative error allowed for each 
dependent variables. The CSP-derived (A2.2) can be numerically integrated using some non-stiff 
solver using the CSP-recommended integration step-size. At. The error that can be tolerated by g 
per integration step is therefore: 8g ■ Ay err or/At. Hence the theoretically neglected g fast is 
numerically negligible during the integration process whenever lg fas N < 5g, which serves as the 
criterion for declaring gfast exhausted. When gfast is declared exhausted, (A2.3) yields. 


f m = B k {F^-F*} =0, m = 1, 2, ... M, 


(A3.1) 


which are the M equations of state referred to earlier. In general, f™ - 0 is achieved by near 
cancellation of the positive and negative terms. Certain terms participate strongly in this balancing 
act, while others are more or less irrelevant. The participation index, P k , is introduced to measure 


the degree of participation: 


pttt 

F ±k“ 


R m F (k) 

B k F ± 


IB^ffv+F^I + I b m •• 8g I 


(A3. 2) 


where the k in numerator is not summed and the r in the denominator is summed over the 173 
reactions. As defined, the maximum magnitude of P™ k is 0.50. In the present study, terms in 
(A3.1) are ordered in descending magnitude of its P^ . and sufficient terms are kept so that the 
the positive and negative sums are at least 0.48 and -0.48, respectively. 


Similarly, not all reactions in (A2.4a) contribute equally to g slow . The importance index, 
l! , is introduced to assess the degree of importance of the k-th reaction to the i-th element of y : 
c i F (k) 

4 . = ^ , (A3.3) 

^ tcX+F^ + lSg 1 ! 


the k and r indices on the right hand side are treated as before. Again, terms in (A2.4a) are order in 
descending magnitude of 4k s0 that the sum of I ily, I is above some user-specified value (the 
present paper used 0.9). 
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In addition, CSP associates with every exhausted mode one or more unknowns which are 

the "fast" variables responsible for the rapid decay of the mode. The diagonal dimensionless 

elements of the NxN matrix a b^ m (no summation on m), denoted by Q(m; i), are called the 

m 

radical pointers of the m-th exhausted mode on the i-th species. The i-th species is said to be a fast 
variable or a radical for the m-th exhausted mode i/Q( m; i) is not a small number in comparison to 
unity. The species with the largest Q(m; i) is identified as the radical for the m-th exhausted mode 
in this paper. If duplication occurs, the next largest pointer is used. The M relations ^(y) ® 0 can 
be used to solve algebraically for these M radicals in terms of the others. If the "wrong species" 
are identified as radicals, the solutions of the M relations are not accurate 14 and can not be trusted. 
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Figures 

Fig. la. Mass Fractions of CH3OH and CO2 vs. time. 
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Fig. lb. Mass Fractions of CO and H vs. time - showing comparison of results of full and 

minimum sets discussed in §7. 
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Abstract 

When very fast chemical reactions are present in the study of re- 
acting flows, the so-called partial-equilibrium and the quasi-steady 
approximations are frequently applied to simplify the chemical kinet- 
ics source terms, to remove stiffness, and to reduce the number of 
dependent variables. An important consequence of these approxima- 
tions is that the mass diffusion term of the individual species must 
be simultaneously modified. This paper presents a general and sys- 
tematic derivation of the modifications, including that of the bound- 
ary conditions, for massively complex chemically reacting flows using 
computational singular perturbation (CSP). Simple examples shall be 
used to demonstrate the theoretical concepts. 


1 Introduction 

In reacting flow problems, when vastly disparate chemical reaction time s- 
cales are present, the partial-equilibrium and/or quasi-steady approximations 
are frequently employedM’®’®’^’® to simplify the chemical reaction source 

•This work is supported by NASA Langley’s Aerothermodynamics Branch, Space Sys- 
tems Division. 
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term. Theoretically, the mass diffusion terms in the PDE’s are concurrently 
affected by this simplification®. In the present paper, we shall present a 
unified derivation of the simplified equations, paying particular attention to 
the modified mass diffusion terms, including the modified boundary and ini- 
tial conditions, under the assumption that the characteristic diffusion Peclet 
number of the problem is moderately large — i.e. the flow is primarily con- 
vection dominated, with diffusion being important mainly in relatively thin 
boundary layers. 

The present work is an extension of the theory of computational singu- 
lar perturbation (CSP)^’®’®’^. While the existing CSP theory has been 
limited to a system of ODE’s, the present extension includes spatial diffu- 
sion and deals for the first time with a system of PDE’s. Anticipating that 
research in modern reacting gasdynamics will be strongly CFD oriented, 
some special numerical issues will also be addressed. It will be shown that 
the CSP methodology of refinement of basis vectors is directly applicable to 
this class of PDE problems, allowing the algorithmic treatment of massively 
complex convection-reaction-diffusion problems which contain exhausted fast 
reactions. 

We shall first work through a simple example using conventional method- 
ology, then present a formal extension of the CSP theory. Finally, the CSP 
methodology is applied to the example as a demonstration. 


2 The General Statement of the Problem 


We consider a reacting gasdynamics problem with N chemistry related un- 
knowns. Let y be the iV-dimensional column vector representing the un- 
knowns. The dot product operator of this iV-dimensional space shall be 
denoted by ”0” to distinguish it from the ordinary of the 3- dimensional 
spatial dot operator. The general dimensional PDE system for yean be ex- 
pressed as follows: 


Dy 

Dt 


g + d, 


( 1 ) 


where g(y) is the chemical reactions source term, d is the mass diffusion 
term: 


d = V* (D © Vy), 


( 2 ) 
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and D is a N x N matrix of diffusion coefficients (usually diagonal). Appro- 
priate initial and boundary conditions for all N unknowns are provided. 

The chemical reaction source term g is usually expressed as the sum of 
R terms each representing an elementary reaction: 

g = f>f\ (3) 

r = l 

where s r and F r are the stoichiometric vector and the reaction rate for the 
r-th elementary reaction, respectively. In general, N ^ R, and in most 
practical situations R is significantly larger than N. However, it should be 
obvious that only N stoichiometric vectors are linearly independent. Hence, 
the right hand side of (3) can always be expressed in terms of N linearly 
independent s r ’s serving as the “basis vectors” of the JV-dimensional space. 

In the modern age of CFD, solutions to such non-linear PDE problems 
can be obtained numerically. The main numerical issues are: the much larg- 
er number of unknowns in comparison to non-reacting gasdynamics, and the 
“stiffness” of the equations introduced by the disparate chemical time scales. 
So long as an adequate computing budget is available, these issues can be 
resolved. However, the amount of physical insights which can be extracted 
from numerical solutions is usually limited. In the present paper, we shall fo- 
cus on the derivation of an approximate simplified non-stiff system of PDE’s, 
including its boundary conditions. The rationale is that the knowledge of 
the approximate simplified equations, in addition to the solution itself, can 
provide much more physical insights than otherwise available. 

2.1 Preliminary Discussion 

When the problem under investigation is sufficiently simple, the convention- 
al approach is to first non-dimensionalize the dependent and independent 
variables, obtain the dimensionaless governing equations, and identify the 
relevent dimensionless Peclet numbers and Damkohler numbers. The initial 
task of non-dimensionalization requires considerable skill; order of magni- 
tude estimates must be made intelligently, and intuition and experience of 
the investigator play a major role. In contrast, the subsequent task is rela- 
tively straightforward: develop an asymptotic theory in the limit of certain 
Damkohler numbers being infinitely large — assuming that all other dimen- 
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sionless parameters remain finite. In most applications, only u leading order 
results are needed. 

In using CFD to compute for the desired solutions, the flow domain of 
interest must first be covered by a suitable mesh or grid system. Let U be 
the local characteristic velocity of the flow field, and h the local characteristic 
cell width. The local characteristic flow time, Tfl ow , is: 

_ h_ 

T flow = jj • 

The local characteristic diffusion time , r, is: 

_ h 2 
T diff= J) 

where D is a representative local diffusion coefficient. The ratio of to 
Tflow t ^ ie l° ca ^ Peclet number: 

m) = ^ ( 6 ) 

which is /i-dependent. In the present paper, we shall assume the Pi{h m ax) 
is moderately large so that globally diffusion does not dominate convection. 

When a chemistry source term g is present, a number of chemical reac- 
tion time scales are introduced into the problem. It is obviously that their 
presence should have a significant impact on the grid system used for the 
computations. At the present time, no general theory exists to guide the 
selection of the appropriate local values of h. Usually, the grid is chosen by 
the investigator based by some informed expectations of the behaviors of the 
desired solutions. As a general principle, one would like to use as few grid 
points as possible without sacrificing accuracy and spatial resolution. Usu- 
ally, a global upper limit h max is placed on h so that the resulting numerical 
solution is assured of adequate spatial resolution. 

The issue at hand is how to deal with reacting flow problems in which 
the characteristic time scales of some fast reactions are very, very small. 


( 4 ) 

( 5 ) 
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3 A Simple Example 

Consider a simple reaction system consisting of two elementary reactions and 
two reactants A and B: 


Reaction #1: 
Reaction: #2: 


A + A «-► B, 
A ^ B. 


(7a) 

(7b) 


The mass diffusion coefficients of A and B are denoted by D A and Db, 
respectively. The dimensional governing PDE for the two reactants are: 


DA 


Dt 

DB 


Dt 


-2 F 1 - F 2 + V- (D a VA), 
F 1 + F 2 + V- (D b VB), 


(8a) 

(8b) 


where D/Dt 


is the substantial derivative operator: 


D 

Dt 


d_ 

dt 


+ v- V, 


(8c) 


v is the flow velocity, and 

F 1 = ki(A 2 - K X B), (8d) 

F 2 = k 2 (A — K 2 B). (8e) 

For the sake of simplicity, we shall assume that the rate coefficients k\, k 2 , 
equilibrium constants K\,K 2 are constants, the flow is steady, and that the 
characteristic Peclet number of the problem is moderately large so that (8a) 
and (8b) are essentially parabolic. At the upstream boundary, initial condi- 
tion is provided. At the "side” and downstream boundary surfaces, appro- 
priate boundary conditions are provided. 

We are interested in deriving the simplified PDE’s and their modified 
boundary conditions for this simple problem when one of the chemical reac- 
tions is very fast. 

For the sake of concreteness, we shall assume that reaction #1 is the fast 
reaction. 

The chemical reaction system also introduces time scales. When confront- 
ed by a large and complex chemical reaction system, however, the above the- 
oretical task is simply not viable. Instead, theoreticians usually proceed by 
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intelligently applying ad hoc approximations which must be guided by intu- 
ition and experience. The so-called partial equilibrium approximation (PEA) 
requires that the theoreticians can somehow identify the fast reactions, the 
so-called quasi-steady approximation (QSSA) requires that the theoreticians 
can somehow identify the fast reactants (frequently called “radicals”). These 
conventional techniques will be demonstrated in the next sections, with spe- 
cial emphasis on the diffusion terms. 

If M (linearly independent) reactions are assumed exhausted, the net re- 
action rates of these M reactions are set to zero. In addition, M species are 
judiciously identified as “radicals” and they are solve for in terms of the re- 
maining N-M unknowns from the M algebraic equations obtained by PEA. 
The subsequent derivation of the leading order simplified equations requires 
some care, but is relatively straightforward. If M species can somehow be 
identified as radicals, the so-called quasi-steady state approximation (QSSA) 
is available. To apply QSSA, the convection term of the radicals are simply 
negelcted. 

If diffusive effects are absent, the simplification of this system of ODE’s 
can be handled by exisiting methods. The focus of this paper is to show how 
diffusive effects should be handled. 

For this simple problem, the dependent and the independent variables 
can be intelligently non-dimensionalized, and the appropriate asymptotic 
theory can be formally developed, the dimensional parameters are: L, the 
characteristic length scale, U, the characteristic velocity of the flow field, 
and D a , D b , h, h and K i, while K 2 is dimensionless. The characteristic 
convection time scale tjj ow is: 

T flow = UV. ( 9 ) 

We shall assume that all the diffusion coefficients have comparable magni- 
tudes so that they can be represented by a single characteristic value D. The 
characteristic diffusion time T d iff * s f^ en: 

T diff= L2 ! D - ( 10 ) 

In addition, we denote the characteristic time scales of reactions $1 and $2 
by Ti and Tj, respectively. We shall be interested in the case when 

n « r 2 « Tfl ow < r di ff. 


( 11 ) 



Appendix III 


87 


The Damkohler number is T di ff/r chern . Then the Peclet number P e for the 
problem, defined as the ratio of Tdijf to tji ow , is: 


Pi = 


UL 
D ' 


( 12 ) 


By assuming that Pi is 0(1) or larger, we exclude the case where diffusion 
may be the fastest process in main domain of the problem. The issue at 
hand is to derive a simplified governing PDE for this example problem when 
reaction $1, which is somehow known to be much faster than reaction $2, 
is already “exhausted.” The modifications to the diffusion terms and the 
modified boundary conditions are of particular interest. 


3.1 The Partial-Equilibrium Approximation 

The conventional derivation of the simplified PDE’s for the example problem 
using PEA would proceed as follows. 

Reaction #1 is somehow identified as the fast reaction, and is assumed to 
be exhausted after a brief transient period (away from boundary surfaces). 
The exhaustion of reaction #1 implies that the net reaction rate F 1 is small 
(but not zero) compared to either the forward or reverse rates which must 
therefore balance each other approximately. By setting 

F 1 = ki(A 2 - KiB) « 0, (13a) 

which is an equation of state, an algebraic relation between A and B. We 

can solve B in terms of A , or vice versa: 

« 4 -, ( i3b ) 

« y/KvB. (13c) 

It is important to emphasize that, in the PEA procedure, (13a) (and therefore 
(13b) and (13c)) must never be substituted into (8a) or (8b)— because (13a) 
is in fact obtained from (8a) or (8b) by taking the formal limit of k x -► oo. 
See (17) later. 

Adding (8a) and two times (8b), we obtain, without approximation: 

^-(A + 2 B) = F 2 + V- ( D a VA ) + 2V- (D b VB), 


B 

or 

A 


( 14 ) 
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which does not contain F l . Goddard [6] called A + 2B the reaction invariant 
of reaction # 1. 

We can now eliminate either B or A by using either (13b) or (13c), re- 
spectively. Using (13b) in (14) to eliminate B, we obtain: 


DA K x 


-F 2 + 


I< x 


Dt K\ + 4A K\ + 4A 
where Dab is given by: 


V- (D ab VA), 


Dab = Da 4- 4A-^— , 


and F 2 is approximately 


.2 ~ hEi A (- x 
~ *» ( /<2 




Alternatively, using (13c) in (14) to eliminate A, we obtain: 


DB 


2 VB 


F 2 + 


4 y/B 


Dt ~ y/TTx+ty/B yfK[ + \\[B 


jV- (D ba VB), 


(15a) 


(15b) 

(15c) 


(16a) 


where is given by: 

^B.4 = Db + (16b) 

and F 2 is approximately 

F 2 « k 2 K 2 y/B(^- - VB). (16c) 

A 2 

Equations (15a) and (13b), or (16a) and (13c) which are their equivalents, 
are the simplified PDE’s, obtained with the PEA method, valid in the flow 
field where the “infinitely fast” reaction #1 is exhausted. They are not valid 
in thin layers immediately adjacent to boundary surfaces where reaction #1 
may be active. Hence, the initial and boundary conditions for (15a) or (16a) 
must be modified. In particular, they are required to satisfy the equation of 
state, (13a). 

By comparing (8a) and (15a), we find: 


F 1 ~ 


K\ + 2A 2 

Kr+4A 


ifrrn v -<^-]ra v - (DsVB) (17) 
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which is also consistent with (8b) and (16a). Note that (17) does not contra- 
dict (13a). Mathematically, it is a more accurate version of (13a) — the terms 
on its right hand side are the next “higher order” terms which were neglected 
in (13a) in comparison to the terms which contribute to F 1 on the left hand 
side. Equation (17) shows clearly that when reaction #1 is exhausted, its 
net reaction rate F 1 is not ”zero;” in fact, its leading order value can be 
accurately expressed in terms of the currectly active chemistry and diffusion 
terms. 


3.2 The Quasi-Steady State Approximation 

The conventional derivation of the simplified PDE’s for the example problem 
using QSSA would proceed as follows. Either A or B must first be identified 
as a “radical” somehow, and then the convective term in the corresponding 
PDE is neglected in comparison to the chemical source and the diffusion 
terms. 

3.2.1 The QSSA applied to B 

Neglecting the convective term in (8b), we solve for F 1 to yield: 

F 1 « -F 2 - V- ( D b VB ). (18) 

Unlike (13a), Equation (18) is not an algebraic equation relating A and B\ 
instead, it is a PDE for B. Using (18) to eliminate F 1 from (8a), we obtain: 

— « F 2 + V- (D a VA) + 2V* (D b VB). (19) 

Dt 

In order to make further progress, we must take advantage of the fact that 
reaction #1 is known to be fast: the forward and reverse reaction rates 
contributing to F 1 on the left hand side of (18) are much larger than the 
terms on its right hand side. Thus, (18) is consistent with (13a), and therefore 
(13b), to “leading order.” Using (13b) to eliminate B from (19), we obtain: 

n A 

^ « F 2 + V- (D ab VA) (20) 

where F 2 is given by (15c). Comparing (18) and (20) with (17) and (15a), 
it is seen that QSSA when applied to B agrees with PEA only when the 
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additional assumption 

Ki » 4A (21) 

is satisfied. 

3.2.2 The QSSA applied to A 

Neglecting the convective term in (8a), we can solve for F 1 to yield: 

F'*J-F 2 + ±V-(D a VA). ( 22 ) 

Using similar procedures as above, we obtain: 

^ « \f 2 + V- (D ba VB) (23) 

where F 2 is given by (16c). Comparing (22) and (23) with (17) and (16a), 
it is seen that QSSA when applied to A agrees with PEA only when the 
additional assumption 

K X «4A (24) 

is satisfied. 

3.2.3 Comments on The Conventional Derivations 

The above analyses can formally be justified when k\K\ » k 2 I< 2 , A = 
0 (Ki/K 2 ) and B = 0(Ki/K$), and that the characteristic Peclet numbers 
are 0(1) or larger. When the order of magnitude estimates are satisfied, the 
PEA is valid when K 2 = 0(1), the QSSA applied to A is valid when K 2 << 1, 
and the QSSA applied to B is valid when K 2 U 1. If the order of magnitude 
estimates are not satisfied, then the validity of the results becomes uncertain 
even when the inequality k\K\ » k 2 K 2 is guaranteed. 

For massively complex realistic problems, it is not possible perform the 
formal asymptotic analysis. In its stead, the rationales justifying the approx- 
imations used must come from the accumulated intuition and experience of 
the investigators. While the validity of the approximations used can in princi- 
ple be assessed by estimating or calculating the next higher order corrections, 
this is seldom done in practice. 
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3.3 The Modified upstream Condition 

If the boundary condition specified at the upstream boundary is not con- 
sistent with the exhaustion of reaction #1, a thin “initial” layer will exist 
to allow the fast chemistry effects to run its course while the diffusive ef- 
fects are expected to be negligible. The standard CSP “radical correction” 
procedure^! can be used to obtain the modified initial condition near the 
upstream boundary. 

3.4 The Modified Side Boundary Conditions 

The simplified PDE’s derived above are valid (after the brief initial layer) 
over the whole physical space of interest except for thin layers adjacent to 
boundary surfaces. The actual boundary conditions specified on the bound- 
ary surfaces must be modified to yield the proper boundary conditions to be 
applied at the edge of the thin layers. 

Let T) denote the physical coordinate normal to the boundary surface, and 
£ denote the physical coordinates which lie on the boundary surface. Within 
a thin layer adjacent to a boundary surface, the PEA and QSSA are not 
expected to be valid. Intuitively, the substantial time derivative term and 
the contribution of reaction #2 to the chemistry source term are expected 
to be small here — in comparison to that of reaction # 1 — and are therefore 
negligible. In other words, reaction #1 and diffusion balance each other 
here. Let the “leading order approximation” to A and B in this thin layer 


be denoted by Am and B b i. We have: 

o * + (25a) 

0 - +F ’ + |< P ^)’ (25b) 

which are now ODE’s. The boundary conditions are: 

Aw (£,»7 = (M) = A 0 ((,t), (given), (26a) 

B b i(Z,r) = 0,t) = B 0 ((,t), (given), (26b) 

A 2 bl ((,oo,t) « K x B hl (t, oo,t), (26c) 


where (26c) is applied at the edge of the thin layer symbolically located at 
1 ] — ► oo. 
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The "boundary” value of A near the boundary surface (77 = 0 + ) is iden- 
tified with the value of A bl at the edge of the thin layer (t? = 00 ). In other 
words, the boundary condition for (15a) is given by: 

A((,rj ~ 0 + 5 *) ** A/(£>o°»0* ( 27 ) 

In the general case when D A and D B depend on A bl and B bl , numerical 
solution of (25a) and (25b) is needed to obtain the values at the edge of the 
thin layer. 

For the special case when D A and Db are constants, explicit analytical 
results can be obtained. Adding (25a) and two times (25b), we have: 

0*-^(D A A bl + 2D B B bl ), (28) 

and its analytic solution which remains bounded as r} — » 00 is: 

D A A b i + 2D B B b i = D a A 0 + 2 D b B q . (29) 

Evaluating (29) at the edge of the thin layer and using (26c), we obtain 
A b i(C,oo,t): 

A b i(£,oo,t) « — + 8 ~Db{D a A 0 + 2DbBo)JKi — D A ]. (30) 

4 u B 

which depends on A 0 , B 0 , D A and D B - The corresponding value of B b i(£, 00 , t ) 
is given by (26c). 

4 Extending CSP to Include Diffusion 

For a reacting flow problem with a massively complex chemistry source term, 
progress can only be made if the investigator has sufficient experience and 
insights to identify the appropriate PEA and QSSA approximations. Even 
after such identification is made, The above derivation is viable only for 
relatively simple systems because of the amount of algebra involved. 

We shall extend the existing CSP methodology, which so far can only 
deal only with a system of ODE’s, to convection-reaction-diffusion problems 
which are governed by a system of PDE’s. The goal is to provide an unified 
algorithm to derive the simplified PDE’s for massively complex problems. 
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4.1 The Basic Idea of CSP 

The basic idea of CSP is to divide the iV-dimensional g-space into a fast and 
a slow subspace. 

4.2 The Fast and Slow Subspace 

Following the existing CSP theory, the column vector g, the chemical source 
term, is formally divided into a fast and a slow components: 

g = g iati + g* ,ow . (31) 

The fast component g fa3i resides in a M-dimensional fast subspace which is 
spanned by a set of trial (column and row) basis vectors a m and b m satisfying 
the orthonormal condition: 

b m 0a n = C m,n = 1,2, , M. (32) 

The question of how the trial fast basis vectors are to be chosen will be 
addressed in the next sections. The value of M is determined by the user- 
specified desired spatial resolution of the solution. For example, if the de- 
tailed structure of the solution of length scale below A is not of interest, then 
M is determined by the requirement that the local velocity times the slowest 
of the fast time scales in g fa,t be less than A. In reacting CFD codes, A 
would be the local grid size. 

Using the fast subspace projection matrix , Q(M), defined by: 

M 

Q (M) = £ a m b m , (33) 

m=l 

we can express g* att in terms of M linearly independent modes: 

g fa3t = Q(M)Og (34) 

M 

= £ a «r. < 35 > 

m=l 

where f m is the amplitude of the m-th mode: 

f m 2 b m © g, m = 1,2, , M. 


( 36 ) 
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The slow component g 3,ow is simply the complementary component of g: 

g .iow = (i _ Q(M)) © g. (37) 

The question to be answered is: what happens to f m (and therefore to g^ a,< ) 
when the fast reactions are “exhausted?” 


4.3 The Asymptotic Solution for f m 

To answer the above question, we need to know what happens to f m following 
a parcel of fluid. Taking the substantial derivative of / m , we have: 1 

r) fm M 

tyr= E A n (/" - /») + b m O J © d > m = l,2,...,M, (38) 

D* m=1 

where 

J = ^ , (N x N Jacobian matrix), (39a) 

dy 

Dh m 

A™ = (— + b m 0 J)0a n , m,n = 1,2,..., M, (39b) 

fZ = -b"© (I-Q(M)) 0 g, m,n = 1,2,..., M. (39c) 

Note that both the M x M matrix A” and the amplitudes /™ depend on 
the fast basis vectors chosen. It is clearly highly desirable for to be as 
“small” as possible. 

The theory of CSP generates from any reasonable set of trial basis vec- 
tors a new set of refined basis vectors with the desirable consequences. The 
refined (row) fast basis vectors, b™, are given by the CSP step #1 refinement 

procedure:$bUO] 

M nun 

K = E'Ct-nT + b"© J )> m = 1,2,..., M, (40a) 

n=l 

and t™ is the inverse of A™. The refined (column) fast basis vectors, a^, are 
given by the CSP step #2 refinement procedure: 

M n~ 

< = E(- i KT + 3 0«.K», rn=l,2,...,M, (40b) 

n= 1 Ul 

1 The considerable amount of algebra involved in the derivation of (38), (39c) and (40a) 
is straightforward and is omitted here. 
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where r£ m is the inverse of A" m . 

Note that these refined basis vectors auto- 

matically satisfy: 




br© a„ = 

m , n — 1,2, . . 

-,M, 

(41a) 

and 

Ko K = C. 

m,n = 1,2,. 

. . ,M. 

(41b) 


In essence, the refinement procedure described above is simply an extension 
of the so-called power method 1111 for finding the largest (left and right) eigen- 
vectors of a real matrix by iteration. 

We shall use the CSP convention that superscript and subscript o indicate 
that the appropriate CSP-refined basis vectors are used®. For example, the 
new Q °(M) is given by: 

M 

<K(AO = E <W- («) 

m=l 

The eigen-values of the M x M matrix t™ have the dimension of time, 
and can be interpreted as the time scales of the fast modes. In the theory of 
CSP, all M fast time scales are assumed asymptotically small in comparison 
to the currently active time scale. Let r(M) denote the slowest of the fast 
time scales, and let r(M + 1) denote the fastest of the slow time scales which 
is of the same order as the fluid mechanics time scale. The parameter e(M) 
defined by: 

e(M) = 

is a measure of the time scale separation between the fast and slow subspace. 
The present theory, following the CSP theory for ODE’s, formally assumes 
large time scale separation, t.e. e(M) << 1. 

The right hand side of (1), g + d, can now be decomposed by using the 
refined basis vectors: 



Dy 

Dt 


= g +*= +<£'"■)+ 


+ d 


o,slow 

o 


), 


(44) 


where 


g o, fast + d o, fast = Q°(M) © (g + d), 


(45a) 
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M 

= E <C(/7 + br O d ). («•>) 

m=l 

P = b™0g, m = 1,2, ... ,M, (45c) 

g o,*/ou; + d o,./ow = (I - Q°(Af)) © (g + d), (45d) 


The issue at hand is to find the asymptotic solution for (g°’ /o4t + d°' /ost ), and 
to show that it can be neglected in (44) after a brief initial transient. 

Using (31), (35) and (37) in (45c), the fast mode amplitudes p can be 
expressed a s: 

/; = r-c ( 46 ) 

Using (46) to eliminate f m from (38), we obtain: 


DP 

Dt 


M 


E A ” (y? - (Co - ■>:© d »- 


m=l 


(47) 


where 


P 

J 0,00 


£ T m( 


m=l 


DJZ 

Dt 


+ 


D b m 
Dt 


O d). 


(48) 


In contrast to /™, which is not guaranteed to be small in the small t(M) 
limit, fo^oo * s formally 0(e(A/)) in the small e(Af) limit provided that d is of 
order unity. 

Assuming that all eigen-values of A™ are essentially negative, we can 
conclude that the leading order asymptotic solution to (47) is: 


P » Co - *C© d, m = l,2,...,M, (49) 

which is valid after a brief initial transient (following the trajectory of a fluid 
parcel). Using (49) in (44) and (45b), we obtain: 


and 


Dy 

Dt 


M 


E <CC + ( J - <KM) © (8 + d >- 


m=l 


(50a) 


ojaat i jo, fast 
o o ' 


M 

E <C 


m= 1 


(50b) 


Both (50b) and (50a) are valid after the brief initial transient layer. 


<?' 7 ~ 
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By neglecting in the small e(M) limit from (50b) and (50a), we 
obtain: 

^ « (I-Q:(M))0(g + d). (51a) 

gO,/a»t + jo./oit ^ 0, (51b) 

The solution of (51a) automatically satisfies (51b) for t>t 0 in the Lagrangian 
sense, provided that it is satisfied at t — t 0 . 

By neglecting f™ >00 in the small e(M) limit from (49), we obtain: 

/- ~ _b-Q d, m = 1,2, . . . ,Af. (52) 


However, except inside a very thin reaction-diffusion layer adjacent to bound- 
aries, (b"0 d) in (52) is itself a “higher order” term in comparison to the 
individual terms which contribute to /“. Hence external to this thin lay- 
er, (b”0 d) can be neglected to obtain M algebraic relations between the 
components of y: 

/ 0 m « 0, m = 1,2,..., M, (53) 


which are called equations of state in the theory of CSP. However, this higher 
order term (b£0 d) must be kept in (52) in the derivation of (51b) and (51a) 


above. _ 

A variant of (50a) was derived by Maas and Pope 1121 who called their 

approach the method of intrinsic low-dimensional manifold (ILDM)^]. In- 
stead of the CSP-refined basis vectors, they used the local fast eigen- vectors 
of J to evaluate the fast subspace projection matrix. In addition, they chose 
not to use (37) for the evaluation of g°’ t,ow , and handled the right hand side 
of (50a) differently. A discussion of their approach can be found in published 
comments which followed the paper in reference [12]. 

The remaining question to be resolved is the modified boundary condition 

for (50a). 


4.4 The Reaction-Diffusion Layer 

Immediately adjacent to physical boundaries, the user-specified boundary 
conditions may not be consistent with (53), and the fast modes are in gen- 
eral active there. A very thin layer therefore exists in which mass diffusion 
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balances fast chemistry. We shall denote the solution in this thin layer by 
y bh the local coordinates on the boundary surface by £ and rj as defined 
previously. 

In this reaction-diffusion layer, we rewrite (1) as follows: 


d = -g + 


Dyu 

Dt ' 


The diffusion term d can be approximated by: 

d - |< D 0 If >• 


(54) 


(55) 


The right hand side of (54) can be decomposed into its fast and slow com- 
ponents using the refined basis vectors: 


S ID «& 


(-g +%):•'““+ <-g + 


Dt 


Dt 


(56) 


Neglecting the slow component and the convective term in the fast compo- 
nent, we obtain the leading order governing equation for the thin reaction- 
diffusion layer: 

£(D 0 ^) » (5T) 

Or) dr) 

Equation (57) is a system of N second order ODE’s and it poses a two- 
point boundary-value problem. The “initial” value for y w at i) = 0 is known 
from the original boundary condition. Since (57) is invariant under the trans- 
formation t) — > —T), “decaying” and “explosive” modes will appear in pairs. 
The CSP initial-value algorithm can be adapted to handle this special class 
of “infinite horizon” boundary layer problems! 14 ! — the explosive modes can 
be explicitly suppressed using CSP methodology. 

The “thickness” of the layer, Atj, is of the order of 0(^D max T(M)) where 
D max is the largest eigen-value of D. It can be considered a “thin” layer 
whenever Arj « L where L is the characteristic length scale of the given 
physical problem. The edge of the layer, denoted by t) — ► oo, is located 
at the point when all M fast modes are declared exhausted. The solution 
of the thin layer equation, (57), remains constant beyond this point, and is 
identified as y ;,/(£, oo, t). 
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The boundary condition for (50a) is: 

y(C v = o + > 0 = yw(£> °°> *)■ ( 58 ) 

It is perhaps interesting to note that the value of y w (£,oo,<) cannot be ob- 
tained by the so-called radical correction using the refined basis vectors 
obtained for the PDE’s — unless all mass diffusion coefficients are identical. 

4.5 The Initial Transient Layer 

In the initial transient layer, the effect of mass diffusion is generally negligible 
and therefore can be treated by the radical correction as demonstrated in 
reference [10]. 


5 The Simple Example Revisited 

To get started, we can use either the fast eigen-vectors of J or the chemical 
stoichiometric vectors of the (guessed) fast reactions as our trial fast basis 
vectors. For our simple example, we consider reaction #1 to be the fast 
reaction: 


a, = 1-2, if, b 1 = [-1,-1]. 
The CSP step #1 refinement yields: 


b 1 = 

° K x + 4 A 


[-2 A,K X \. 


The CSP step #2 refinement yields: 


i, =s [-2 + 0(e(l)), 1 + 0(£(l))f. 


(59) 


(60) 


(61) 


where the time scale separation parameter c(l) is proportional to k 2 /k x (K x + 
4A). It is seen that the refined a? can be approximated by the trial ai in 
the limit of k x — ► oo when e(l) is a small parameter. This will always be the 
case whenever the trial column basis vectors have been “judiciously” chosen. 
In fact, the closeness of a m and a^ can be used in general as a “test” of how 
good the original choice of the trial vectors were. In practice, one should 
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proceed with the CSP analysis only when a m is a good approximation to a^. 
It can easily be verified that if the trial basis vectors were randomly chosen, 
the refined basis vectors will rapidly converge — if c(M) is small — after a few 
full cycles of CSP refinements. 

For the sake of simplicity, we shall omit here the 0(e(l)) term in the fast 
subspace projection matrix: 


q:(i) 


i 


K,+4A 


4A 
-2 A 



(62) 


The corresponding slow subspace projection matrix is: 


i-q:(i) 


i 


Ki + 4A 


( 


K i 
2A 


2/Ci \ 

4 A ) * 


(63) 


The analytical developments in §3 which assumed reaction #1 is fast can 
readily be reproduced by appropriately using these two projection matrices. 
For example: (15a) is the first component of (51a), (13b) is (53) with M — 1, 
(25a) and (25b) are “leading order” versions of (57) with N — 2. 


6 Discussion 

Mathematically, the above general theoretical development is an asymptotic 
analysis of a singularly perturbed problem in the limit of e(M) — ► 0 where 
e( M) is the ratio of the slowest exhausted fast time scale to the current time 
scale. Unlike the original CSP theory which is capable of including higher 
order terms by recursive applications, the present extension for PDE’s is 
limited to only the leading order approximation. In return for this concession, 
the governing ODE’s in the very thin reaction-diffusion layer are cleanly 
decoupled from the PDE’s governing the rest of the flow field. 

The value of /” as given by (48) should be monitored in the course of 
the computation. Its value should be “small,” and may be included in (50a) 
and (50b) to improve the accuracy of the solution. If this small correction 
term is found to be significant, then one or more of the exhausted fast modes 
is reviving and the value of M should be decreased. In general, the updating 
of M along a fluid streamline can be dealt with in the standard manner. 
When an energy equation is included in the original model, then the term 
representing thermal conduction must also be appropriately modified^. 
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In a viscous boundary layer near a solid boundary, the flow field velocity 
is small as a consequence of the no-slip condition, and increases rapidly with 
rj. In addition, the temperature profile in a thermal boundary layer is often 
also rapidly varying. As a consequence, the value of M and the identifica- 
tion of the fast reactions are expected to vary significantly across the viscous 
and thermal boundary layers. In the very thin reaction-diffusion layer, how- 
ever, the flow velocity does not enter and the temperature is expected to 
be approximately the wall temperature. While theoretically the identifica- 
tion of the fast reactions can change across this layer, no strong variation 
is expected. In the present presentation, we have tacitly made the ad hoc 
assumption that the reduced reaction system representing at the edge 

is valid throughout the very thin reaction-diffusion layer. 

An important assumption in the present theoretical development is that 
all eigen- values of A™ are essentially negative. In other words, the fast modes 
are assumed to be decaying modes. If A™ does contain positive eigen-values, 
then the use of (49) will also suppress the fast “explosive modes.” Such 
modes are physically essential and must be “unleashed” to describe certain 
important processes such as the ignition of flames. In addition, when A” 
contains essentially imaginary eigen-values, then the use of (49) will also 
suppress the fast “oscillatory modes.” The proper treatment of fast explosive 
or oscillatory modes which must not be suppressed is a difficult matter and 
is not yet adequately addressed by CSP at this time. 
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Abstract 

The Computational Singular Perturbation (CSP) method of sim- 
plified kinetics modeling is reviewed with emphasis on the justification 
of the methodology. A simple example is first worked through using 
the conventional methodology of partial-equilibrium and quasi-steady 
approximations, and then treated in some detail using CSP. 


1 Introduction 

When one is confronted with an unfamiliar problem in chemical kinetics, the 
traditional first step is to identify the relevant chemical species and the impor- 
tant elementary reactions which occur among them. The resulting “complete 
model” of the reaction system is usually further simplified to take advantage 
of available approximations. For sufficiently simple problems, conventional 
analytical methods can be usedM’^’^W 4 !-^!. In most situations, the success 
of the resulting simplified model is measured not only by its quantitative 
predictive capability, but also by its simplicity — the fewer superfluous terms 
the better. Generally speaking, approximate analytical results are highly val- 
ued because of the insights they can provide when inspected by a competent 

•This work is supported by NASA Langley’s Aerothermodynamics Branch, Space Sys- 
tems Division. 


105 



Appendix IV 


106 


theoretician. But when the reaction system is massively complex, this is not 
a practical option. 

Recently, databases containing extensive, reliable and up-to-date data 
for certain reaction systems are available. Computations using complete 
models from such databases can now be routinely carried out. In this new 
computational era, it is no longer necessary to pick out only the relevant 
chemical species and the important elementary reactions because the pres- 
ence of benign superfluous terms in the formulation is not a problem. An 
option increasingly available to the modern theoretician is to first generate 
a complete model numerical solution, examine the resulting data to discern 
significant and interesting causes and effects, making additional diagnostic 
runs if necessary, and then try to propose simplifications and approxima- 
tions. Why do theoreticians still care about simplified approximate models 
when double-precision numerical solutions to the complete model are easily 
available? The reasons are: physical understanding and “stiffness.” A theo- 
retician would like to be able to make general statements about the problem 
in addition to showing color slides of the numerical solutions. The observed 
behaviors of the computed solutions need to be described in terms of fa- 
miliar concepts, such as chain-branching, chain-termination, ignition delay, 
building up of the radical pool, heat release etc. In particular, theoreticians 
would also like to be able to identify the rate-controlling reactions (for the 
chemical species of interest in the time interval of interest), fast reactions for 
which rate coefficients do not need to be known accurately, and superfluous 
reactants and reactions which need not be included at all. In addition, the 
vast disparity of time scales which is responsible for the simplifications and 
approximations is also responsible for stiffness^, a generally undesirable at- 
tribute of the governing differential equations from the computation point of 

The theory of computational singular (to 

be referred to as CSP here) exploits the power of the computer to do simpli- 
fied kinetics modeling. In essence, CSP is a systematic mathematical proce- 
dure to do boundary- layer type singular perturbation analysis. While it can 
be used to obtain analytical results for simple problems, it is designed to be 
used for massively complex problems using computations. A CSP computa- 
tion not only generates the numerical solution of the given problem, but also 
the simplified equations in terms of the given information. Most interesting 
questions about the reaction systems can be answered merely by inspection 
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of the numerical CSP data. 

The present paper is a review of the theory of CSP. The basic CSP re- 
sults are summarized, and are applied to a simple example to illustrate the 
essential features. 


2 The Mathematical Problem 

We shall consider a spatially homogeneous reaction system consisting of R 
elementary chemical reactions. The total number of unknowns, which in- 
clude concentrations of the chemical species and other state variables such 
as temperature and total pressure, is N. We shall represent the N unknown- 
s by a ^-dimensional column vector, y = [j/ 1 ,?/ 2 , ■ • • ,y N ] T - The governing 
equation for y is a system of quasi-linear ODE: 

ft = (2 ' 1) 

where g is the sum of contributions from the R elementary reactions: 

g = f> r F r (y), (2.2) 


and s r and F r ( y) are the (generalized) stoichiometric vector and the reac- 
tion rate of the r-th elementary reaction, respectively. The value of R may 
be greater, less than, or equal to N. We shall call (2.2) the physical repre- 
sentation of g, because each additive term can be satisfactorily explained by 
the investigator who formulated the problem. 

Usually, an investigator is interested only in certain special species, each 
with a different accuracy requirement and perhaps in a different time in- 
terval. In many cases, not all the initial conditions needed to compute the 
solutions are known, and many rate constants needed in the computations 
are only educated estimates. The mathematical problem is to derive the sim- 
plest model of the reaction system consistent with the user-specified accuracy 
requirements. 
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3 Some Definitions 

When the forward and reverse reaction rates of a single or a group of fast 
reactions are in approximate balance, we say the reaction or the reaction 
group is in partial-equilibrium. When the production and destruction rates 
of a particular species are in approximate balance, we say the species is 
in quasi-steady state. In either case, an approximate algebraic relation is 
obtained between components of the y vector. Such relations, which do not 
contain free parameters, shall be called equations of state. In many cases, 
species in quasi-steady state are chemical radicals, but not always. 

When M equations of state are available, they can be used to selectively 
solve for M unknowns in terms of the others — in so doing the need for M of 
the ODE’s is eliminated. Which M unknowns should be solved for in terms of 
others, and which M ODE’s should be discarded? Lam [12] provides explicit 
theoretical guidance to make these choices, and calls the M unknowns CSP 
radicals, or simply radicals for short, when they can be accurately solved 
for from the AI equations of state. When used in this CSP context, the 
term radical carries a special meaning distinct from the context of chemical 
structure. In most situations, a CSP radical is also a chemical radical, and 
vice versa. But it is not always true. The identification of CSP radicals is 
done via a “radical pointer” which shall be demonstrated later. 

Frequently, additional (exact or approximate) algebraic relations between 
the components of the y vector may exist. The conservation law of atomic 
species is one such example. Such relations are distinguished from the above 
equations of state by the presence of free parameters which are determined 
by initial conditions. Following Lam [12] , we shall call all such relations con- 
servation laws. 


4 A Simple Example 

We shall use a simple hypothetical reaction system^ 2 ^ to illustrate the issues 
involved. 

Let the state vector be y = [A, B, C] T where A and B are chemical 
concentrations and C is temperature. The reaction system consists of three 
elementary reactions: 

reaction# 1 : A + A ^ B, 


(4.1a) 
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reaction # 2 : A *=* B, (4-lb) 

reaction #3 : B + B ^ A. (4.1c) 

The (generalized) stoichiometric vectors and the reaction rates are: 

Sl = [-2, 1 , A H,} T , F 1 = ki(A 2 - I<iB), (4.2a) 

s 2 = [-1, 1, A H 2 ] T , F 2 = k 2 (A - K 2 B), (4.2b) 

s 3 = [1, -2, A H 3 ] T , F 3 = k 3 {B 2 - K 3 A). (4.2c) 


where the reaction rate coefficients k\,k 2 ,k 3 and the equilibrium constants 
k^k 2 ,k 3 are known and — for the sake of simplicity — their dependence on 
C is assumed negligible. The heats of reaction (with the proper units) for the 
three reactions are denoted by A H u A H 2 and A H 3 , respectively. We shall 
find it useful to separately identify the forward and reverse reaction rates as 
follows: 


F r = F r + -FL, r = 1,2, 

(4.3) 

where F+ and FL are both positive. 
The system of ODE is: 


-jr — S\F l + s 2 F 2 + s 3 F 3 , 

at 

(4.4) 

which can be written out as follows: 


^ = —2 F 1 - F 2 + F 3 , 
dt 

(4.5a) 

= F 1 + F 2 — 2F 3 , 
at 

(4.5b) 

— = AH.F 1 + AH 2 F 2 + AH 3 F 3 . 

(4.5c) 


To make things concrete, the rate coefficients are given numerical values: 


k\ fa 

1 0 4 cc/ mole- secon d , 

I < , 

« 1.1 x 10 2 mole/cc, 

(4.6a) 

k 2 ~ 

10 _1 /second, 

I < 2 

« 1.1 x 10 2 , 

(4.6b) 

k 3 « 

10 4 cc/mole- second, 

k 3 

fa 0.8 x 10~ 8 mole/cc, 

(4.6c) 


and 


A Hi 

fa +1.1 x 10 4 cc-°F/mole, 

(4.7a) 

A H 2 

fa +1.0 x 10 5 cc-°/<'/mole, 

(4.7b) 

A H 3 

ra —2.9 x 10 5 cc-°F/mole. 

(4.7c) 
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The initial conditions are also given numerical values: 

>1(0) « 1.5 x 10 -4 mole/cc, 5(0) « 0.1 x 10 -6 mole/cc, C(0) « 300° if. 

(4.8) 

The investigator is interested in A(t), accurate to two significant figures, after 
the first few seconds up to a few minutes. 

Experience and intuition can play no role here because the problem is 
hypothetical, and indeed may not even make “chemical sense.” Note that 
detailed balance would require K 3 = I<i/(K 2 ) 3 , thermodynamics would re- 
quire A H 3 = AHi — 3AH 2 , and the Law of Mass Action would require s r 
and F T to be consistent. The data provided do not satisfy some of these re- 
lations exactly — they were chosen to demonstrate that the techniques under 
discussion are completely mathematical, and do not require consistency with 
constraints “external” to the given initial-value problem. 

Because this example is sufficiently simple, a conventional asymptotic 
analysis can be performed — provided that the dependent and independent 
variables can be intelligently non-dimensionalized, and a small dimensionless 
parameter can be identified. The reader can readily confirm that even for 
this simple problem the task of non-dimensionalization of variables is not 
straightforward. Consequently, the conventional method usually proceeds 
with a dimensional formulation without the benefits of dimensionless param- 
eters. Progress is made via intuitive ad hoc judgment on the speed rankings of 
the elementary reactions, and the use of partial-equilibrium or quasi-steady 
approximations must follow certain special procedures. We shall illustrate 
the conventional methodology in the following section. 


5 The Conventional Methodology 

We shall demonstrate the partial-equilibrium approximation first, followed 
by the quasi-steady approximation. 

5.1 The Partial-Equilibrium Approximation 

For the example problem, a competent investigator will (correctly) conclude 
that reaction #1 is fastest, reaction #2 is next, and reaction #3 is slowest. 
In the time period of interest, it is expected that reaction #1 will have 
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exhausted itself, reaction #2 is the rate- controlling reaction, and reaction 
#3 is essentially dormant. 

Reaction #1 is exhausted: The partial-equilibrium approximation for re- 
action #1 consists of setting the net reaction rate of reaction #1 to zero. 
The following equation of state is obtained: 

F 1 = h(A 2 - I<iB) 0, (5.1) 

Note that (5.1) contains no free parameters. One may use it to solve 
A in terms of B, or B in terms of A. 

Equation (5.1) must be handled with care — it must never be substi- 
tuted directly into the original system of equations. If this advice is 
not heeded and (5.1) is substituted into (4.5a), (4.5b) and (4.5c), the 
resulting equations are simply wrong. This is because the small net 
reaction rate of an exhausted fast reaction is in general competitive 
with the currently active slower reactions. In other words, (5.1) is only 
adequate to be used as an equation of state, but is inadequate to be 
used in the original equations. See (5.10) later. 

Reaction #3 is dormant: The contribution by reaction #3 to (4.4) is neg- 
ligible and need not be included in the time period of interest. It is 
emphasized that reaction $3 being negligible does not mean F 3 ~ 0. 
Hence, B 2 « K 3 A is not valid and must never be used. By neglecting 
s 3 F 3 from (4.4), one can show that: 

(A Hi - A H 2 )A + ( AHi - 2A H 2 )B + Constant. (5.2) 

which is a (temporary) conservation law. An astute investigator will 
be able to attribute the physical meaning of “conservation of total 
energy” to (5.2). If one were interested in the reaction system over a 
time period of the order of hours, reaction #3 would not be negligible 
and (5.2) would not be valid. 

The Simplified Model: The procedures recommended by Williams^! to 
follow-up the partial-equilibrium approximation proceeds as follows. 
First, one of the ODE in the system of equations is used to eliminate 
F 1 from the rest of the equations. For example, (4.5b) can be used to 
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eliminate F 1 from (4.5a) and (4.5c). These two equations are then sup- 
plemented by (5.1) differentiated with respect to time. The simplified 
model is obtained by solving for dA/dt, dB/di and dC/dt from these 
three equations. We obtain: 


dA 

„ K ' 

3^1 F 3 

(5.3a) 

dt 

~ K l+ iA F 

A", +4/1 ’ 


dB 

_ 2A F i 

M f* 

(5.3b) 

dt 

~ K,+4A F 

K\ + 4A ’ 


dC 

dt 

» [AH 2 — Ai/i ( 

Kl +2y4 'llF2 

Ki+4A )]F 



+ [A^3 + A/ fl ( 2 ^ + + ^ ) )]F3. 

(5.3c) 


Note that either (5.3a) or (5.3b) can be replaced by the algebraic e- 
quation of state, (5.1). Since F 1 does not appear in this simplified 
model, this system of ODE is no longer stiff (regardless of how large 
jfcj is). It can easily be verified that its solution satisfies (5.2) when 
( AHi - 3A H 2 - A H 3 )F 3 is negligibly small. 

The New Initial Conditions: The initial conditions at t — * 0 + are no 
longer given by (4.8), since they must satisfy (5.1). A detailed analysis 
will show that reaction #2 can also be considered dormant in the brief 
initial transient period. Hence, in addition to (5.2), A + 2 B is also 
approximately conserved is this period. Using these results, we have: 

A(0 + ) « A(0) + 2(5(0) -5(0+)) = 1.46 x 10" 4 , (5.4a) 

5(0 + ) « j 94 x iQ- 6 . (5.4b) 

K\ 

The value C(0 + ) « 301.04 can be found from (5.2) applied between 
t = 0 and t = 0 + . 

5.2 The Quasi-Steady Approximation 

Alternatively, an investigator may choose to proceed with the quasi-steady 

approximation instead of the partial-equilibrium approximation. 
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The quasi-steady state approximation requires that certain species be 
chosen as “radicals.” If B is chosen to be the radical, we neglect dB/dt from 
(4.5b) to yield the following equation of state: 

F 1 » - F 2 + 2 F 3 . (5.5) 


Unlike (5.1), this equation is substituted directly into (4.5a) and (4.5c) to 
eliminate F 1 . We obtain: 

^ « F 2 — 3F 3 , (5.6a) 

dt 

^ « [AH 2 - AHJF 2 
at 

+ [AH 3 + 2AH 1 ]F 3 , (5.6b) 

These equations agree with (5.3a) and (5.3c) but only in the K\ » A limit. 
The radical, B , is to be solved from (5.5) in terms of A. 

UKt « A, then A should be chosen as the radical instead, yielding: 

F 1 « (-F 2 + F 3 )/ 2. (5.7) 


Substituting this equation into (4.5b) and (4.5c) to eliminate F 1 , we obtain: 
(4.5c) to eliminate F 1 . We obtain: 


dB_ 

dt 

dC 

dt 


^ F 2 ^ F 3 
2 F ~2 F ' 

[A H 2 - iAff.lF 2 
+{AH, + iA«,]F\ 


(5.8a) 


(5.8b) 


These equations agree with (5.3a) and (5.3c) but only in the K\ « A limit. 
The radical, A, is to be solved from (5.7) in terms of B. 

Whenever K\ — 0(A) , the quasi-steady approximation simply cannot be 

use df 12 M 13 l. 

According to this presentation, the quasi-steady approximation is seen to 
be a special case of the partial-equilibrium approximation. See §5.4 later. 
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5.3 How Good Are These Results? 

The analytical results obtained above are only the leading order term of an 
asymptotic theory valid in the limit of k\ — ► oo. How good are they for the 
“real” problem at hand? For the example problem, the small dimensionless 
parameter c for an asymptotic theory turns out to be either k^Ki/^iKi 
or k 3 /(kiK%), whichever is larger. For the given rate data, we have c = 
kiK-ifk^Kx « 10-\ which is barely adequate to provide one significant figure 
in the answers. More seriously, if the order of magnitudes of the initial data 
>1(0) and B( 0) are significantly different from (4.8), these asymptotic results 
may be totally misleading. 


5.4 The Exhausted Fast Reaction Rate 


We have shown that the equation of state obtained by the partial-equilibrium 
approximation, (5.1), must not be substituted directly into the original equa- 
tions because it has inadequate accuracy, while the equation of state obtained 
by the quasi-steady approximation, (5.5) or (5.7), can be more freely used 
with the caveat that its validity depends on the relative magnitude of A 
versus K\. The question is: can a better approximation be found for F 1 for 
large time than either (5.1), (5.5) or (5.7)? 

Differentiating F 1 with respect to time, we obtain: 

, /ft dA __ dB . 

k ' (2A n - *•-*> 

-h(Ki -I- AA)[F l - F^], (5.9a) 


dF 1 

dt 


where 


Ki+2A 2 2 (K x +A) f3 

Ki+tA K^ + tA 


(5.9b) 


It is now clear what F 1 does as time marches on. Initially, F 1 decays expo- 
nentially with characteristic time scale 0(|r(l)|), but eventually it follows F^ 
which evolves with a slower time scale. The long time asymptotic solution 
of (5.9a) is, for t » 0(|r(l)|): 


F 1 


+ r(l) Jt Fl + 


1 

ki(Ki + 4 A) 


(5.10) 


where 


t(1) = 


(5.11) 
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In other words, the leading approximation for F 1 is F when reaction #1 
is exhausted. If F 1 « is substituted directly into (4.5a) and (4.5c), we 
will indeed recover (5.3a) and (5.3c), the partial-equilibrium results obtained 
earlier. Hence, F 1 « F^ is valid without restriction on the relative magni- 
tudes of K\ and A\ it is thus a superior equation of state than either (5.1), 
(5.5) or (5.7), and agrees with them in the appropriate limits. 

6 The Theory of CSP 

The present section explains the CSP procedures and provides intuitive jus- 
tifications for their validity. 

6.1 Observations on the Conventional Methodology 

The most important information needed by the conventional methodology is 
the speed ranking of the reactions and the identification of the CSP radical- 
s. Once the fast reactions are somehow identified, and the radicals chosen, 
the partial-equilibrium and the quasi-steady approximations are available to 
make further progress — with appropriate caution on the unreliability of the 
latter. The derivation presented in §5.4 is a new and different way to make 
progress — it does not need to identify the radical and it clearly yields the 
most accurate equation of state, F 1 F^. All it needs is assurance that 

reaction #1 is the fastest. 

The theory of CSP provides a systematic way to identify the fast reac- 
tions. In addition, it generalizes the procedures in §5.4 to find the exhausted 
feist reaction rates for massively complex reaction systems. 

6.2 Basis Vectors 

The vector g contains all the physics of the problem, and is usually given by 
the investigator formulating the problem using the physical representation, 

( 2 - 2 )- . . 

In general, an N-dimensional vector may be expressed in terms of any set 
of N linearly independent basis vectors^!. CSP exploits the theoretician’s 
prerogative to express g in an alternative representation , and look for basis 
vectors with special properties. 
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Let a i(t), i — 1, 2, . . . , N, be a set of N linearly independent column basis 
vectors. The set of inverse row basis vectors, b'(t), i = 1,2 ,...,N, can be 
computed from the orthonormal relations: 

b' 0a i = Sj, ij = 1,2,...,N, (6.1) 

where © is the dot product operator of the ^-dimensional vector space. 
The column vector g can now be expressed as: 

g = £>/', (6-2) 

t=l 


/* = b* © g = £ BtF r , > = 1,2 JV, (6.3) 

r=l 

and 

B* r = b* © s r , i = 1,2, . . . ,N. (6.4) 

Each of the additive terms in (6.2) represents a reaction mode, or simply 
mode. The amplitude and direction of the i-th mode are /* and a;, respec- 
tively. Eventually, CSP provides an algorithm to find the “best” set of basis 
vectors for the derivation of the simplified models. 

The physical representation of g uses the physically meaningful (and time- 
independent) stoichiometric vectors as the default column basis vectors. For 
our example, the default set is: 


ai = Si = [—2, 1, AH\] t , (6.5a) 

a 2 = s 2 = [-l, 1, A H 2 ) t , (6.5b) 

a 3 = s 3 = [1, -2, A H 3 ] T . (6.5c) 

Using this set, the inverse row vectors can easily be computed: 

b 1 = [2AH 2 + AH 3 , AH 2 + AH 3 , 1 }/H, (6.6a) 

b 2 = [—2AHi — AH 3 , —AH\ — 2AH 2 , -3 }/H, (6.6b) 

b 3 = [ -AHx + AH 2 , — AH\ 2AH 2 , -1 ]/H, (6.6c) 

where 

H = AHi — ZAH 2 — AH 3 . (6.6d) 



Appendix IV 


117 


Using the given input numerical data, we have: 

H = 10 3 cc-°A'/mole. (6.7) 

Inspite of this quite respectable dimensional value, H is actually very nearly 
“zero,” because it is a small number in comparison to AHi, AH 2 and AH 3 . 
For “sufficiently small” H, an alternative a 3 must be provided. We shall take 
advantage of H ± 0 here and proceed with caution. 

It is straightforward to verify that at t = 0, we have: 



P = 
p = 
P = 

b 1 0 g = F 1 = 2.14 x 10 -4 mole/cc-second, 
b 2 0 g = F 2 = 1.39 x 10 -5 mole/cc-second, 
b 3 © g = F 3 = — 1.19 x 10 -8 mole/cc-second. 

(6.8a) 

(6.8b) 

(6.8c) 

6.3 

The Speed Ranking of the Modes 


Differentiating (6.3) with respect to time along a solution trajectory y (t), we 
obtain: 

C=EA i = (6.9) 

dt 

where 

a; 

= f + b-'Qj] © aj, i,j = 1,2,... ,N, 

at 

(6.10a) 


j 

Qgr 

= ^ N x N Jacobian matrix. 

dy 

(6.10b) 


The non-linear nature of the original problem is manifested by the fact that J 
is in general not a constant matrix. At any moment in time, the eigen- values 
of J can be computed. For problems arising from chemical kinetics, they are 
usually essentially real — when the problems are of the boundary-layer type. 
The reciprocal of an eigen-value, called the time scale, has the dimension of 
time, and shall be denoted by r(i). Ordering them in increasing magnitudes, 
we have: 

|r(l)| < ... < |r(t)| < ... < |r(iV)|, (6.11) 

which provides an approximate speed ranking of the “eigen- modes.” 

The question what is an ideal set of basis vectors now has an obvious 
answer: ideal basis vectors should diagonalize Aj, thereby uncoupling all the 
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modes. This is the standard strategy for analyzing linear problems, but it 
obviously needs some modifications for non-linear problems because there is 
no finite algorithm to diagonalize A}. A full discussion of this subtle point will 
be given in a future paper. When the set of basis vectors used is non-ideal, 
the modes are coupled, and each mode will not have a distinct characteristic 
time scale. As shown in our example in §5.4, the fastest mode F 1 evolves 
with its own characteristic time scale |r(l)| only initially. As it becomes 
exhausted, it eventually follows F^, which evolves with the characteristic 
time scales of the slower modes. This mode mixing is the price we pay for 
not diagonalizing A j, and is an intrinsic issue of non-linearity which must 
be dealt with: the fast modes behave as slow modes when they are near 
exhaustion. Viewed in this light, the task of deriving simplified models is 
reduced to finding basis vectors such that the fast modes mix with the slow 
modes as little as possible. From the pragmatic point of view, however, 
it is not really necessary to uncouple all the modes; it will be sufficient to 
uncouple the fast modes from the slow modes approximately — so that the 
residual coupling can be neglected in accordance with some user-specified 
accuracy threshold. Mode mixing among the fast modes or among the slow 
modes causes no difficulty and can be tolerated. 

6.4 The Classification of Fast and Slow Modes 

First of all, we need a more precise classification of fast and slow modes. 
Usually, an investigator has a definite idea on the desired time resolution At 
of the solutions — events whose time scales are below At are not of interest. 
Hence, the group of M modes which satisfy: 

|r(m)| < At, m = 1,2,...,M, (6.12) 

are considered fast modes , and all others are the slow modes. The fastest 
group of active slow modes are the rate- controlling modes. Slow modes with 
negligible amplitudes are called dormant modes. 

If one is interested in all time scales, then At = t(M + 1)- When time is 
measured in units of At, the fast modes should all be nearly exhausted. We 
shall deal with a precise definition of exhaustion later. In §5.4, we showed how 
to obtain the “asymptotic” large time solution for F 1 . We shall generalize 
the procedure to an arbitrary fast mode below. 



Appendix TV 


119 


6.5 The CSP Refinement Strategy 

CSP does not attempt to diagonalize A’ . It assumes that at any moment 
in time an intelligent set of trial basis vectors is available, and provides an 
algorithm to refine it. The strategy is to provide a systematic, programmable 
algorithm which generates a new A’ which is more block-diagonal than before. 

For non-linear problems, the eigen-vectors of J are time-dependent, and 
thus they do not diagonalize Aj. They are, however, excellent trial basis 
vectors. As demonstrated in the previous sections, the conventional method 
uses the default stoichiometric vectors and requires a good guess of the speed 
ranking of the modes. 

The CSP theory uses the ratio 


= 


r( M) 
t(A/ + 1) 


(6.13) 


as a small dimensionless parameter, and develops the refinement algorithm 
by an asymptotic analysis in the small cm limit. Physically, is a measure 
of the time scale separation of the fast and slow modes. Each application of 
the CSP refinement procedure will depress the magnitude of the off-diagonal 
blocks of A j by O(cm)- 


6.6 The CSP Refinement Procedure 

We shall assume that at any moment in time on a solution trajectory, the 
value of M is known, and a set of trial basis vectors for the fast modes are 
given: 

a m , b m , m — 1,2, ... ,M. (6-14) 

With only the fast basis vectors available, we can compute for the M x M 
upper-left block of A* , denoted by u>™: 


,db 




= A™ = (— — + b m © J) © a n , m,n = 1,2,... ,M. 


dt 


The inverse of shall be denoted by r™: 


M 

M 


= E t >: 

n'= 1 

n'=l 


(6.15) 


m, n = 1,2, . . . , M. 


(6.16) 
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We shall mark the refined basis vectors and all entities evaluated with refined 
basis vectors by either a superscript or a subscript o, or both. 

To refine the b m row vectors, we use the so-called step # 1 procedure: 

M db n 

K = EC(j + b n 0J), m = l,2,...,M, (6.17a) 

= a m , m = 1,2,...,A/. (6.17b) 

The step #1 procedure depresses the upper-right M x (N — M) block of 
A’ — it makes the fast modes “purer” by weakening their couplings with the 
slow modes. To refine the a m column vectors, we use the so-called step #2 
procedure:: 

b™ = b m , m = 1,2, ... ,M, (6.18a) 

M da 

a° m = E(-^+ J m = 1,2,. ..,M. (6.18b) 

n=l 

The step #2 procedure depresses the lower-left (N — M) x M block of A J- — it 
makes the slow modes purer. Note that at each refinement step the orthonor- 
mality condition is always satisfied. These two steps may be performed singly, 
or in tandem in any order, or recursively any number of times — provided the 
most current t™ is used always. This refinement procedure is mathematical- 
ly equivalent to that presented in Lam and Goussis^, the difference being 
that in this formulation the slow basis vectors are not involved. If the time 
derivative terms in (6.17a) and (6.18b) were omitted, the procedure would be 
identical to the so-called Mises Power Method for finding eigen-vectors asso- 
ciated with the largest eigen-values (See Carnahan, Luther and Wilkes^!). 
In essence, by allowing mode mixing among the fast modes, CSP extends the 
Mises Power Method to compute the next iterant for the fast basis vectors. 
Since the time scales of the fast modes are all faster than the current time 
scale, the time derivative terms are always small corrections. 

In practice, the first guess of the trial basis vectors is usually time- 
independent. The subsequent refined sets will in general be time-dependent 
because they are, by construction, y-dependent. In any case, their time 
derivatives can be evaluated accordingly. In a computer program, these time 
derivatives can be evaluated approximately using either stored or predicted 
data already available in the integration routine. Such programming issues, 
however, are beyond the scope of this paper. 
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6.7 The Fast Subspace Projection Matrix 

We can form a N x N matrix Q(M) as follows: 

M 

Q(M) = J2 a - bm > (6-19) 

m=l 

and call it the fast subspace projection matrix. When evaluated with refined 
basis vectors, it will be denoted by Q£(M). 

Any column vector u or row vector v can be decomposed into its fast and 


slow components using either Q (M), Q 0 (M) or Q°(M). 

Using Q°(M), we 

obtain: 

u = u° 0 ' fast {M) + uy l0W (M), 

(6.20a) 


v = v° 0 ' fast (M) + w°/ ow {M\ 

(6.20b) 

where 

u°' fa3t (M) = Q ° 0 (M) 0u, 

(6.21a) 


u 0 o ’ ahw {M) = (I-Q° 0 (M)) ©u, 

(6.21b) 


v°/ a3i (M) = v © Q° 0 (M), 

(6.21c) 


v o,si° W ( M ) = v © (I — Q°(M)). 

(6.21d) 


We can decompose Q(M) into its M components: 

M 

Q (M) = £ Q m , (6-22) 

m = l 

where 

Q m = a m b m , m = 1,2,... ,M. (6.23) 

We shall call Q m the fast mode projection matrix of the m-th mode. 

The radical pointer of the m-th mode, Q m (i), is given by^^: 

Q m (i) = the i-th diagonal element of Q m , i = 1, 2, . . . , N. (6.24) 

Note that Q m (i ) is dimensionless, and its sum over all N components is 1.0. 
Geometrically, Q m (k ) is a measure of the projection of the k-th unit vector 
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in the m-th mode. Hence, whenever Q m (k ) is not a small number, species k 
is said to be a CSP radical. 

The fast reaction pointer of the m-th mode, P m (r), is given by: 

Pm(r) = (s r ) _1 O Q m O S r , r = 1,2,...,*, (6.25) 

which is a dimensionless number nominally of order unity. Here, (s r ) _1 is a 
set of row vectors orthonormal to the s r set of (linearly independent) column 
vectors. Physically, P m {r) is a measure of the projection of the r-th stoichio- 
metric vector in the m-th mode. Whenever P m (r) is not a small number, the 
r-th reaction is said to be a fast reaction. 


6.8 The Main CSP Results 

Using the available refined basis vectors, the original system of ODE’s be- 
comes: 

% = g °S ait (M) + g°/ ow (M), (6.26) 


dt 


where 


n 


= q :(M) eg = EK:'‘" F ' 


R 




r=l 


r=l 


M 


(6.27a) 

(6.27b) 


K:l‘“ = Q:(M) 0 s P = E *C B Z ■ r=l,2,...,«, (6.27c) 


m = 1 


R 


/ 7 = K 0g = EV 171 = 1,2, ... ,A/, 

r=l 

B™ r = b;0s r , m— 1,2,. . . ,Af, 


and 


g r low (M) = (i-qs(m)) 0g = Es:f”-r, 


(6.27d) 

(6.27e) 

(6.28a) 


r= 1 


s° o y™ = (I - Q SCAT)) G s r = s r - £ < 5 o m r7 (6.28b) 
r = 1,2, ... ,& 


M 

E 

m=l 
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We shall call s 0 o ,a J 0w and s° 0 \{ a,t the effective stoichiometric vectors of the r-th 
reaction. 

The unrefined fast mode amplitude f m satisfies: 

ifm M 

% r = '£»Z(r-&), m = l,2,...,M, (6.29) 

n=l 

where 

fZ = (b m -b-) 0g = / m -/; (6.30a) 

= -b;oErt m = 1,2, , M, (6.30b) 

r=l 

and in (6.29) and in (6.30b) are evaluated using the trial basis 
vectors. Equation (6.29) is the generalized version of (5.9a). 

Treating (6.29) as a linear equation for / m , we can express its solution as 
the sum of a homogeneous solution and a particular solution: 


r = + fZ’t- < 6 - 31 ) 

We shall assume that the eigen- values of are essentially negative (i.e. 
the fast modes are of the boundary-layer type). The homogeneous solution 
is responsible for satisfying the initial conditions, and for boundary-layer 
type problems it becomes exponentially small in time measured in units of 
r(M). The total adjustment of y in this brief time period can be given 
approximately by: 

M 

(Ay). c = - £ [w:(n)(r-rj M1 («- 32 ) 


derived under the ad hoc assumption that a m r^ 1 (A/) is approximately time- 
independent. Equation (6.32) is called the radical correction by Lam^ 2] . For 
large time ( t » |r(M)|), f m -> f£ r t.- The asymptotic solution for /£ rt . 
can be shown to be: 


M df n 

= fZ + E C (A*K“ 


n= 1 


dt 


m = 1,2, ... , M. 


(6.33) 
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If g /a,t evaluated with f m ~ f™ ri . is considered negligible in comparison 
to g slow according to some user-specified criterion, the following simplified 
model is obtained: 

^ * g° 0 ’ thw (M). (6.34) 

The “initial” conditions for (6.34) at t = 0 + must satisfy the following ap- 
proximate equations of state: 

fT = K 0g = EV«O, m = l,2,...,M. (6.35) 

r=l 

The radical correction algorithm can be used to enforce this condition as 
shall be demonstrated later. Note that (6.34) provides N ODE’s for the N 
unknowns, and its solutions are guaranteed to satisfy (6.35) for all t > 0+ 
automatically provided (6.35) is satisfied initially and all fast eigen- values of 
J are essentially negative. 

7 The CSP Method on the Example 

We shall step through the example numerically using CSP. 

7.1 Choosing the Trial Basis Vectors 

We assume that at the beginning we have no idea which reaction is fast. The 
eigen- values A(t) and eigen- vectors of J at t = 0 can be computed numerically. 
We have: 

A(l) = -1.27 x 10 2 /second, (7.1a) 

A(2) = -0.173/second, (7.1b) 

A(3) = 0.00/second, (7.1c) 

indicating that there is a fast mode with time scale of the order of 10 -2 
seconds, followed by a slower mode with time scale of the order of about 10 1 
seconds. The right (column) eigen- vectors a, and left (row) eigen- vectors /3* , 
ranked in order of decreasing speed, are: 


<*1 = 

[-1.90, 

0.993, 

1.87 x 10 4 ] r , 

(7.2a) 

a 2 = 

[-0.960, 

-2.46 

x 10" 2 , -0.901 x 10 5 ] r , 

(7.2b) 

a 3 = 

[0.000, 

0 . 000 , 

O 

o 

(7.2c) 
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and 

/3 1 = [-2.46 x 10 -2 , 0.960, 0.00], (7.3a) 

/3 2 = [-0.993, -1.90, 0.000], (7.3b) 

/ 3 3 = [-0.890 x 10 5 , -1.89 x 10 5 , 1.00]. (7.3c) 

These may be used as our time- independent trial basis vectors for t > 0 but 
they diagonalize A .*• only at t = 0. Since our time resolution of interest is in 
seconds, only the first mode can be considered feist. Hence, M — 1. 

Alternatively, we can use the stoichiometric vectors of the three reactions 
as our trial basis vectors (i.e. the default set in §6.2). Since they are time- 
independent, A’- is easily computed: 

/ -1.16 x 10 2 -1.13 x 10 2 2.23 x 10 2 \ 

A‘ = -1.12 x 10 1 -1.11 x 10 1 2.21 x 10 1 . (7.4) 

\ 0.216 x 10" 2 0.208 x 10“ 2 -0.408 x 10" 2 ) 

This matrix has significant off-diagonal terms for all t. Nevertheless, its 
diagonal elements can be used to estimate the time scales of the modes, and 
they indicate that mode #1 (i.e. reaction #1) is fastest and its time scale is 
possibly below the time resolution of interest — suggesting that M = 1. 

In what follows, we shall take M = 1 and choose (6.5a) and (6.6a) as our 
trial fast basis vectors for t > 0: 

ai = [-2.00, 1.00, 1.10 x 10 4 ] r , (7.5a) 

b 1 = [-9.00 x 10\ -1.90 x 10 2 , 1.00 x 10" 3 ]. (7.5b) 

Using these trial basis vectors, the N elements of the radical pointer of mode 
#1 are (taking advantage of the fact that H ^ 0): 

Q x ( 1) = 1.80 x 10 2 , Qi(2) = -1.90 x 10 2 , ^(3) = -1.10 x 10 1 . (7.6) 

Radical pointers computed using trial basis vectors are theoretically unre- 
liable, and this result should be ignored. We shall compute a theoretically 
reliable radical pointer for mode $1 using the refined basis vectors later. The 
R elements of the fast reaction pointer of mode #1 are, as expected: 

P x (l) = 1.00, P x (2) = 0.00, Pi( 3) = 0.00. (7.7) 

confirming that reaction $1 is the fast reaction of the trial fast basis vectors. 
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7.2 Skipping the Rapid Transient Period 

Since we are not interested in the rapid transient period which lasts tens of 
milli-seconds, the main issue now is to find the adjusted initial conditions for 
the simplified model which governs the slow evolutionary period. 

In the rapid transient period, y adjusts rapidly in such a way that the 
amplitude of the fastest mode approaches zero. This adjustment can be 
computed approximately by the radical correction given by (6.32). For the 
example here, the amplitude of the fastest mode at t = 0 is Z 1 = F = 
2.14 x 10 -4 . Making the radical correction using the trial fast basis vectors, 
we obtain the following adjusted initial condition at t = 0 + : 

y(0 + ) = y(0) + Ay rc = [1.46 x 10" 4 , 1.95 x 10" 6 , 300.02] T , (7.8) 

which yields a much smaller amplitude, f 1 = F 1 = —1.36 x 10 -7 . Note that 
the relative magnitude of the correction to B is much larger than those to A 
and C. 


7.3 The CSP Refinement 

Using the trial fast basis vectors (7.5a) and (7.5b), we can decompose g into 
its feist and slow components: 

g = g /<**' 4. g alow , (7.9a) 

where 

g /“ at = Q(l) © g = [2.68 x 10 -6 , —1.34 x 10 -7 , —1.47 x 10 _2 ] r , (7.9b) 
g ,iow _ (I_Q(1)) 0 g = [6.87 x 10" 6 , -6.90 x 10" 6 , -0.667] T . (7.9c) 

In general g faat is of the “same order” as g s,ow (component by component), 
and its neglect will cause an order unity error since it was computed using 
the trial fast basis vector which has never been refined. 

The refined fast basis vectors at t = 0 + are: 

aj = [-1.92, 1.00, 1.88 x 10 4 ] r , (7.10a) 

b* = [-2.52 x 10 -2 , 0.950, 0.000], (7.10b) 

with t°( 1) = (bj 0 J 0 aj) _1 = -0.788 x 10 -2 seconds. Note that a" and 

b 4 are quite close to a x and /3 1 , the right and left eigen- vectors of the fastest 

mode. Note also that b 1 and bj bear no resemblance to each other. 
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Using the above refined basis vectors for t > 0 + , the new A j is: 

/ -1.27 x 10 2 -2.56 x 10" 1 7.70 x 10" 1 \ 

A‘° = 1.03 -1.65 x 10" 1 4.95 xlO" 1 . (7.11) 

\ 3.68 x 10~ 2 9.25 x 10" 4 -2.78 x 10" 3 ) 

Note that the off-diagonal terms are much smaller than before. 

The elements of the new radical pointer using the refined basis vectors of 
mode #1 are: 

Q“,(l) = 0.482 xl0-\ <E,i(2) = 0.952, <£,(3) = 0.00. (7.12) 

This (refined) radical pointer is theoretically reliable, and it indicates that B 
alone is qualified to be a CSP radical, and that C must never be so identified. 
The elements of the new fast reaction pointer are: 

P°i(l) = 0.906, P^(2) = 0.941 x 10"\ J? fl (3) = 3.58 x 1(T 5 , (7.13) 

again confirming that reaction #1 is the fast reaction. 

We can fine tune the initial conditions using the refined basis vectors and 
apply the radical correction once more. We obtain: 

y(0 + ) = [1.46 x 10- 4 mole/cc, 1.89 x 10 _6 mole/cc, 300.02° A'] 7 ’. (7.14) 

This correction is mainly on B , the CSP radical. It is applied again here 
because (7.8) was computed using the unrefined trial basis vectors and was 
only a rough correction. It is possible not to use the radical correction if one 
is willing to numerically integrate the original system of ODE’s including 
g °<f aat and let the rapid transient do the adjustments. 

We can now compute the fast and slow components of g at t = 0 + : 

g°o Jaat = Qo(l) ©g 

= [1.24 x 10" 5 , -6.52 x 10' 6 , -0.128] T , (7.15a) 

g° 0 ' a,ow = (i-q:(i)) ©g 

= [-5.89 x 10 -6 , —1.61 x 10 -7 , — 0.555] T . (7.15b) 

Equation (7.15a) must be interpreted with care — it should not be used to 
assess whether g^ a,t can be neglected in comparison with g* low . The proper 
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order of magnitude estimate of should be made with f™ « fZpart. as 

mentioned earher. The large time asymptotic approximation to the (refined) 
particular solution is: 


M df n 

C„.«C + EC'M J F + -. m=l,2,...,M, (7.16) 

n= 1 


where 


M 


= n =1.2 M, 

n=l 


(7.17) 


and /£ was previously defined in (6.30a). In a computer code, can be 
properly evaluated — the time derivative in (7.16) and (7.17) can be computed 
either exactly or approximately ( e.g . using backward finite difference). A 
rough estimate of the order of magnitude of the exhausted amplitude f™ is: 


ST - fZ 


m 

,oo 


0 { 


Tq{ M ) wm 

r‘(Af + l)" 00 ’ 


m 


1,2,... M, 


(7.18) 


where t°(M ) and r°(M +1) are the time scales of the slowest fast mode and 
the fastest slow mode, respectively. r“(M ) can be estimated by the largest 
reciprocal diagonal element of a A” refined at least once. t°(M + 1) can be 
similarly estimated, or it can be taken to be the integration step size selected 
by an integration routine (e.g. RKQC^^). Hence, we have: 


ojaat 

o o 


M 


M 




’( M ) 


m=l 


m = 1 


r»(A/ + 1) 


)fZ )• 


(7.19) 


For the example, we have t °( 1 ) 0.8 x 10 2 , t°( 2) « 0.6 x 10 1 , and « 

6.63 x 10" 6 . Hence, this rough estimate yields: 

g °/ aat « 0([— 1.6 x 10 -8 , 0.8 x 10" 8 , 1.6 x 10' 4 ] T ). (7.20) 


If all of its components are considered small enough (see next section), we 
can declare the refined fast mode exhausted — inspite of (7.15a) which does 
not appear negligible at all — and neglect g°' faat in comparison to g° 0 '* low to 
yield the desired simplified model. 

For our example, the amplitude of the exhausted mode #1 is: 

11 = b\ 0 g = + b;, 2 F ! + Bl^F 3 . 


(7.21) 
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Neglecting /™ in comparison to the magnitude of the largest positive or 
negative terms on the right hand side, we have: 


EW-P) » 0, (7.22a) 


fX.-F” « 0. (7.22b) 

r=l 

Either is a semi-analytical equation of state — the coefficients B* 0 are provid- 
ed only numerically (at t = 0 + ): 

B l oA = b] © s x = 1.00, (7.23a) 

Bl a = K © s 2 = 0.975, (7.23b) 

b{ z = b l ©s 3 = -1.92. (7.23c) 

The terms on the left-hand side of (7.22a) are nearly balanced, and the dom- 

inant terms can easily be identified by the participation indeJ ^^ — showing 
correctly that Fj_ ~ F\. 

Equation (7.22a) can be used to solve for B , the identified CSP radical, 
thus eliminating the need for its ODE. Equation (7.22b) can be use to solve 
for F 1 , the identified fast reaction, in terms of F 2 and F 3 : 

F 1 = -0.975 F 2 + 1.92F 3 , (7.24) 

which should be compared with (5.1), (5.5) and (5.10), the analytical results 
valid in the limit — » oo obtained earlier. See §8 for additional comments. 

We can now compute for g°' alow : 


„o,slow V' o t slow pr _ o,tloui ril , o,alow pi , o,alow p3 
go — 2L* S o,r r ~ S o,l r ■+" S o,2 r T ®o,2 r > 

r=l 


(7.25) 


where, at t — 0 + , we have: 

8 o,si°w = [_0.916, -0.250 x 10" 1 , -0.863 x 10 5 ] T x 10" 1 , (7.26a) 
8° 0 $° w = [0.859, 0.234 x 10 _1 , 0.809 x 10 5 ] T , (7.26b) 

< 3 °“' = [-2.67, -0.728 x 10" 1 , -2.52 x 10 5 ] T , (7.26c) 
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and 

F 1 = -0.975F 2 + 1.92F 3 , (7.27a) 

F 2 = —6.76 x 10 _6 mole/cc- second, (7.27b) 

F 3 = 2.61 x 10 -8 mole/cc-second. (7.27c) 

Which slow reaction(s) (F 2 or F 3 or both) is controlling g°' ,low can easily be 

identified by the importance index introduced in Lam^l 

The above numerical s° 0 ' a,ow 's derived for the given problem should be 
compared with (5.1), (5.3a) and (5.3c), the analytical results obtained earlier. 
In the ^ -too limit, the CSP-derived analytical results are: 


o t alow 1 

®o,l J k\ —>oo 

[ o t alow i 

S o, 2 Jjbi~»oo 

otalow i 
3 Jfci — *oo 


= [ 0 , 0 , 0 ] T , 
f I<1 2A 
[ Ki+4A' Ki+4A' 
-3 K, -6 A 

h< 1 + AA' + AA' 


A H 2 - A Hi 


Kr + 2A 
Kr+AA 


T 


A H 3 - A ^ 


2{K\ + -d) ,j 
Fj + 4A J 


(7.28a) 

(7.28b) 

(7.28c) 


7.4 The User-specified Accuracy Threshold 

A perceptive reader would have noticed in (7.20) that the error introduced 
to J3, the CSP radical, by the neglect of would be the largest. This 

pessimistic error estimate can be improved by using the radical correction^ 

It can easily be shown that when g° ,faat is neglected and the initial 
conditions properly adjusted, the solution generated automatically satisfies 
f m « 0, m = 1,2,...,M. The theory of CSP uses the radical pointer to 
identify one or more unknowns as CSP radicals for each exhausted fast mod- 
e. The resulting M equations of state can be used to eliminate the ODE’s 
of the M CSP radicals — using the radical correction. Hence, the accuracy 
of the CSP radical is controlled by the equations of state, and not by the 
ODE’s. A detailed discussion of this subtle point is beyond the scope of this 
paper. It suffices to state that the error committed in neglecting g° 0 ' faat can 
be estimated in the conventional way for all other unknowns, but for the 
CSP radicals the correct estimate is a factor + 1) smaller than 

indicated when the radical correction is used. 

With this caveat, the error of neglecting g° 0 Jaat can be computationally 
assessed, and fast modes are declared exhausted only when their neglect in- 
troduces an error estimated to be below the user-specified accuracy threshold , 
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y error . The exhaustion criterion^ ^ for the m-th fast mode is: 

| a m/o m oo T W| « y error, m = 1, 2, . . . , M, (7.29) 

which must be satisfied component by component. 

7.5 Exact or Approximate Conservation Laws 

In a time period of interest, some reaction modes may be so slow that they 
can be considered dormant — they can be neglected. 

Let the (N— M)-dimensional slow subspace be spanned by a j and b* 7 , J = 
M + 1 , . . . , N, the slow trial set which complements the fast trial set of 
basis vectors. The refinement process described earlier will refine these basis 
vectors to weaken the coupling between the fast and slow modes. However, 
within the fast and slow subspaces, the modes are mixed. 

One may look for dormant modes among the N—M slow modes using any 
reasonable algorithm. If the amplitude f* of the K-th (slow) mode satisfies, 
component by component: 

I &0 Kfo T ( M + *)l << y error ) (selected K's), (7.30) 

then it is considered a dormant mode. 

Exact dormant modes occur frequently in chemical kinetics problems — 
they are called conservation laws of atomic species. Generations of students 
have been taught to look for such conservation laws and to use them to 
advantage in analysis. From the CSP point of view, all exact conservation 
laws are simply special cases of dormant modes. However, not all dormant 
modes represent conservation laws. 

Dormant modes play no significant role in CSP — inclusion of dormant 
modes does not cause stiffness. Identification of dormant modes is of interest 
primarily because they may suggest physically interesting concepts. CSP 
provides no special technique to find them; we shall assume that they can be 
somehow found and be identified using (7.30). 

Unlike exhausted fast modes, the relation 

fo= bf © g « 0, (selected K's), (7-31) 

does not produce an equation of state. In the case of conservation law of 
atomic species, for example, b^ © g = 0 is an identity: zero equals to zero. 
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However, if can be expressed in the form: 

bf = (selected /Ts), (7.32) 

where 0 k iO and ® K ,o are scalar functions of y, then mode indeed may 
represent a conservation law. Lam^l called such 0 k>’s eligible conserved 
scalars, and such modes holonomic modes — a name borrowed from classical 
mechanics^ '1. At the present time, there is no known method to deter- 
mine whether a computed slow b* is holonomic — except when it is a time- 
independent constant. All exhausted fast modes are approximately holonom- 
ic, with 0 m , o = 1 and 0 m , o ~ f™. 

Whenever a holonomic /* is dormant, ®k,o becomes a conserved scalar 
which is an algebraic relation among the unknowns. Unlike equations of 
state obtained from exhausted fast modes, conserved scalars always contain 
a free parameter determined by initial conditions. For addition discussion, 
see Lam' 12 !. 

In our example, (7.27c) shows that F 3 is quite small in comparison to 
F 2 . This observation by no means implies that F 3 ph 0 — it only means 
that s °]3 ow F 3 can be neglected from (7.25) in this time period. To find the 
dormant mode, we may use (s x , s 2 , 0 ( 3 ) as our trial basis vectors. The 
slow mode #3 will be found to be a dormant mode, with = [0,0, 1] T and 
b 3 = [-0.890 x 10 s , -1.89 x 10 5 , 1.00]. The row vector b 3 will be found 
numerically to be nearly constant in time, and analytically it can be shown 
that: 


b 3 « [A if, - AH 2 , A H x - 2A H 2 , 1](1 + OihlUlhKJ). (7.33) 

In other words, this is approximately an holonomic mode, with 0 3 , o = 1 and 

0 3 0 = b 3 Oy« (A H\ — AH 2 )A + (A Hi — 2A H 2 )B + C, (7.34) 

which can be interpreted as the total energy of the system. Theoretically, 03, 0 
evolves with time scale of 0(H/(AH 3 k 3 Kz)) which is measured in hundreds 
of seconds; it only appears to be a conserved scalar quantity in the time 
period of 0(l/(k 2 K 2 ) which is measured in seconds. 
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8 A St iff ness- Related Programming Issue 

The Associative Law of Addition is one of the fundamental laws of algebra. 
It states that: 

(x + y) + z = x + (y + z). (8.1) 

It is not generally appreciated that this law is not valid on a finite precision 
computer. One can easily confirm this observation by trying x = 10 15 , y = 
— 10 15 , and z = 1.23 on a hand calculator. 

The numerical violation of the Associative Law of Addition is caused 
by loss of significant figures in subtracting large numbers or adding large 
and small numbers. In a chemical kinetics problem involving fast and slow 
reactions, the evaluation of the vector g indeed involves subtraction and 
addition of large and small numbers. Hence one must pay special attention 
to the evaluation of g inside a computer program. 

Mathematically, g can in principle be evaluated by any of the following 
formally identical expressions: 


R 


g := 

X>F\ 

T— 1 

R R 

(8.2a) 

g : = 

r= 1 r= 1 

M 

(8.2b) 

g := 

where f™ is given by (6.27d): 

£</r+ gr'-". 

r=l 

R 

(8.2c) 

fm 

J 0 

and 

m = 1,2,. . . , Af, 

r= X 
R 

(8.2d) 

g ° 0 ' 3,0W (M) := 

o,»low jpr 
2 1* r * 

r= 1 

(8.2e) 

F r := 

II 

t. 

1 

(8.2f) 


For sufficiently stiff problems, all three evaluations are numerically unreliable. 
However, (8.2c) can be made reliable by applying CSP concepts as shall be 
shown presently. 
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The amount of cancellation in the evaluation of f™ and F r can be mea- 
sured by: 




m 

o 


6F T 


l/i 


m 
o I 


IF'l 


m = 1,2, , M, 


mi+i^r 


r = 1,2, 


(8.3a) 

(8.3b) 


which are both dimensionless numbers. Clearly, when either Sf ™ or 8F T is 
very small, the accuracy of f™ or F r so evaluated is suspect. In our example, 
6f] and 8F 1 are 0(k 2 K2/(kiKi)) a ^ ter t ^ ie ra pid transient period, and — 
to dramatize the situation— if k 2 K 2 l{kiKi) = O(10 -18 ) then even double- 
precision arithmetic would be inadequate. 

CSP uses (8.2c), but intelligently. When /™ is near exhaustion, CSP 
provides its leading order asymptotic approximation, (6.33). Since is a 
known function of y, the needed time derivative can be appropriately eval- 
uated. Hence, whenever f™ has inadequate significant figures as evaluated 
by (8.2d), (6.33) can be used in its place. In fact, it is recommended that 
/ooo used (instead of /”) to assess whether the m-th mode should be de- 
clared exhausted. In other words, CSP properly recognizes the contribution 
°f g°J a3t to g in (8.2c) to be at best a minor correction whenever substantial 
cancellation occurs in the direct evaluation of /™. 

Similarly, in the evaluation of g°' alow , all F r 's are involved. The values 
of the exhausted fast 6F r, s will be very small, and the corresponding values 
of the F r 's will be suspect. The theory of CSP provides M fast reaction 
pointers to identify M fast reactions, and their net reaction rates can be 
solved for from (8.2d) in terms of the rest of the net reaction rates. Hence, 
the use of potentially inaccurate fast F r, s are avoided. In our example, 
the value of F 1 after the rapid transient periods should not be computed 
from F l = fc 1 ( J 4 2 — KiB), the theoretically exact expression, but should 
be evaluated from F 1 « -(^ >2 F 2 + B^F 3 )/ Bl A instead, the CSP-derived 
approximation. 

The computation of the effective stoichiometric vector of the r-th reaction, 
s° 0 '* tow , involves [I - Q°(M)], and its dimensionless diagonal elements are 
obtained by a subtraction process which may lose significant figures. The 
following artifice has been found successful in eliminating specious errors 
in the direct evaluation of s° 0 '* ,ow : Whenever any diagonal element of this 
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matrix falls below an appropriate threshold ( e.g . 10 4 for four significant 
figures calculations), it is replaced by a zero. 


9 Discussion 

The traditional role of a theoretician has always been to simplify seemingly 
complex problems into their bare elements. Physical insights and intuition 
are at the root of this process, helped considerably by the methodology of 
asymptotics when a small parameter is available and can be identified. As a 
practical matter, asymptotic solutions usually consists of very few terms — 
the so-called leading order solution is usually all that could be expected 
because of the massive amount of algebra involved. In §5, the analytical 
results obtained are formally correct only in the limit of ^^ 2/(^1 K\) — ► 
0, kz!(k\Kl) — ► 0. No assurance of accuracy is provided when applied to 
problems with k<iKil{k\K\) « 0.1. 

The theory of CSP welcomes physical insights and experience, but is not 
dependent on them. It uses the eigen-values of J to order the trial modes, 
and provides a refinement procedure to improve the decoupling of the trial 
fast and slow subspaces. From the viewpoint of CSP, users of conventional 
asymptotics simply make educated guesses at the fast basis vectors. The 
special procedures to apply the partial-equilibrium and quasi-steady approx- 
imations are just variants of the CSP refinement procedure. In the absence 
of insights and experience, the conventional methodology cannot get started 
at all, but CSP can proceed routinely — the eigen-vectors of J can always be 
used as the trial set^. Most importantly, the refinement procedure can be 
recursively applied, allowing “higher-order” corrections to be included. The 
CSP user has the responsibility to specify the threshold of tolerable error 
y eTTO r ■, and CSP makes sure that the error introduced by the neglect of 
is below threshold. Different simplified models will be generated for different 
user-specified y error] a lenient threshold will yield a simpler model than a 
more stringent one. 

The method of CSP is a significant advance over the method of matched 
asymptotic expansion for boundary- layer type problems. First of all, it 
can be applied without the need of non-dimensionalizations and identification 
of small parameters. In essence, CSP exploits the disparity of time scales be- 
tween the exhausted fast modes and the currently active slower modes, and 
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the small dimensionless parameter being exploited is e*/ defined in (6.13). 
The simplified model is constructed by an iterative procedure; each itera- 
tion improves the accuracy of the model by 0(e\f)- It is important to note 
that the time derivative term in (6.17a) and (6.18b) must be included in the 
refinement procedures if accuracy beyond “leading order” is desired. The 
superiority of the CSP approach to the conventional methodology, when ap- 
plied analytically to simple problems, is clearly demonstrated in §7.3. For 
massively complex problems, the CSP approach has no peers. 

The CSP-derived simplified model remains a system of N ODE’s which is 
accompanied by a set of M equations of state. Theoretically, the solutions of 
the simplified model automatically satisfy the M equations of state for t > 0+ 
if the initial conditions satisfy them at t = 0 + — assuming that exhausted 
modes did not become active again. CSP uses all N ODE’s to march forward 
in time for all N components of y, and uses the M equations of state only to 
apply the radical corrections to prevent the “drifting” of the exhausted fast 
mode amplitudes. 

For each fast mode, CSP provides a radical pointer which identifies the 
CSP radicals — species which can be solved for from the equation of state. It 
is extremely important to note that one may not arbitrarily select any M 
species to be solved for from the M algebraic equations. In the conventional 
approach, the algebraic difficulty of solving for the radicals frequently forces 
additional ad hoc approximations. CSP deals with this obstacle again using 
iteration; the programmable radical correction procedure can be recursively 
applied to solve for the CSP radicals. In fact, the radical correction should be 
applied after every integration step to counter the usually larger estimated 
error of neglecting g°’7 0it for the CSP radicals. The refined basis vectors 
from the previous time step can be used as trial basis vectors for the current 
time step. In addition, CSP provides a fast reaction pointer for each fast 
mode, and an algorithm to accurately evaluate the amplitudes of exhausted 
fast reactions. 

The CSP method has no difficulty in identifying dormant modes, but does 
not provide a method for finding the so-called eligible conserved scalars — 
except for the simplest case when the row basis vector of the mode in question 
is found to be time-independent. At the present time, CSP does not take 
advantage of any conservation laws which may be available. 

Trevino et. a /.[20],[21] successfully studied ignition phenomena with the 
assistance of CSP data. CSP data can also be used to assess sensitivity of so- 
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lutions with respect to the input rate coefficients^^, and offers an attractive 
alternative to the conventional method of sensitivity analysis^h 

After mode #2 also becomes exhausted, it can be shown by conventional 
methodology that A and B become time-independent, but C varies linearly 
with time when HF 3 ^ 0. The CSP method routinely obtains the simplified 
model in this period with M = 2, and the exercise is left to the interested 
readers. 
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